build diffusion operator for explicit

src_reactDiff/AdvanceDiffusion.cpp

@@ -1,6 +1,9 @@
 #include "reactDiff_functions.H"
 #include "chemistry_functions.H"
 
+#include "AMReX_MLMG.H"
+#include <AMReX_MLABecLaplacian.H>
+
 void AdvanceDiffusion(const MultiFab& n_old,
                       MultiFab& n_new,
                       const MultiFab& ext_src,
@@ -184,5 +187,94 @@ void AdvanceDiffusion(const MultiFab& n_old,
         Abort("AdvanceDiffusion() - invalid reactDiff_diffusion_type");
     }
 
+}
+
+
+void DiffusiveNFluxdiv(MultiFab& n_in,
+                       MultiFab& diff_fluxdiv,
+                       const Geometry& geom,
+                       const Real& time) {
+
+    // fill n ghost cells
+    n_in.FillBoundary(geom.periodicity());
+    MultiFabPhysBC(n_in, geom, 0, nspecies, SPEC_BC_COMP, time);
+
+    BoxArray ba = n_in.boxArray();
+    DistributionMapping dmap = n_in.DistributionMap();
+
+    // don't need to set much here for explicit evaluations
+    LPInfo info;
+
+    // operator of the form (ascalar * acoef - bscalar div bcoef grad) phi
+    MLABecLaplacian mlabec({geom}, {ba}, {dmap}, info);
+    mlabec.setMaxOrder(2);
+
+    // store one component at a time and take L(phi) one component at a time
+    MultiFab phi (ba,dmap,1,1);
+    MultiFab Lphi(ba,dmap,1,0);
+
+    MultiFab acoef(ba,dmap,1,0);
+    std::array< MultiFab, AMREX_SPACEDIM > bcoef;
+    AMREX_D_TERM(bcoef[0].define(convert(ba,nodal_flag_x), dmap, 1, 0);,
+                 bcoef[1].define(convert(ba,nodal_flag_y), dmap, 1, 0);,
+                 bcoef[2].define(convert(ba,nodal_flag_z), dmap, 1, 0););
+
+    // build array of boundary conditions needed by MLABecLaplacian
+    std::array<LinOpBCType, AMREX_SPACEDIM> lo_mlmg_bc;
+    std::array<LinOpBCType, AMREX_SPACEDIM> hi_mlmg_bc;
+
+    for (int idim = 0; idim < AMREX_SPACEDIM; ++idim)
+    {
+        if (bc_spec_lo[idim] == -1 || bc_spec_hi[idim] == -1) {
+            if ( !(bc_spec_lo[idim] == -1 && bc_spec_hi[idim] == -1) ) {
+                Abort("Both bc_spec_lo and bc_spec_hi must be periodic in a given direction if the other one is");
+            }            
+            lo_mlmg_bc[idim] = LinOpBCType::Periodic;            
+            hi_mlmg_bc[idim] = LinOpBCType::Periodic;
+        }
+
+        if (bc_spec_lo[idim] == 0) {
+            lo_mlmg_bc[idim] = LinOpBCType::inhomogNeumann;
+        } else if (bc_spec_lo[idim] == 1) {
+            lo_mlmg_bc[idim] = LinOpBCType::Dirichlet;
+        } else if (bc_spec_lo[idim] != -1) {
+            Abort("Invalid bc_spec_lo");
+        }
+
+        if (bc_spec_hi[idim] == 0) {
+            hi_mlmg_bc[idim] = LinOpBCType::inhomogNeumann;
+        } else if (bc_spec_hi[idim] == 1) {
+            hi_mlmg_bc[idim] = LinOpBCType::Dirichlet;
+        } else if (bc_spec_hi[idim] != -1) {
+            Abort("Invalid bc_spec_hi");
+        }
+    }
+
+    mlabec.setDomainBC(lo_mlmg_bc,hi_mlmg_bc);
+
+    // set acoeff to 0and bcoeff to -1
+    mlabec.setScalars(0., -1.);
+
+    acoef.setVal(0.);
+    mlabec.setACoeffs(0, acoef);
+
+    for (int i=0; i<nspecies; ++i) {
+
+        // copy ith component of n_in into phi, including ghost cells
+        MultiFab::Copy(phi,n_in,i,0,1,1);
+
+        // load D_fick for species i into bcoef
+        for (int d=0; d<AMREX_SPACEDIM; ++d) {
+            bcoef[d].setVal(D_Fick[i]);
+        }
+        mlabec.setBCoeffs(0, amrex::GetArrOfConstPtrs(bcoef));
+
+        MLMG mlmg(mlabec);
+
+        mlmg.apply({&Lphi},{&phi});
+
+        MultiFab::Copy(diff_fluxdiv,Lphi,0,i,1,0);
+
+    }
 
 }

src_reactDiff/reactDiff_functions.H

@@ -26,6 +26,11 @@ void AdvanceDiffusion(const MultiFab& n_old,
                       const Real& time,
                       const Geometry& geom);
 
+void DiffusiveNFluxdiv(MultiFab& n_in,
+                       MultiFab& diff_fluxdiv,
+                       const Geometry& geom,
+                       const Real& time);
+
 ////////////////////////
 // In AdvanceTimestep.cpp
 ////////////////////////

src_reactDiff/reactDiff_namespace.H

@@ -1,73 +1,20 @@
 namespace reactDiff {
 
-    // 0=D + R (first-order splitting)
-    // 1=(1/2)R + D + (1/2)R (Strang option 1)
-    // 2=(1/2)D + R + (1/2)D (Strang option 2)
-    // -1=unsplit forward Euler
-    // -2=unsplit explicit midpoint 
-    // -3=unsplit multinomial diffusion
-    // -4=unsplit implicit midpoint
+    // see reactDiff_functions.cpp for descriptions and default values
     extern AMREX_GPU_MANAGED int temporal_integrator;
-
-    // only used for split schemes (temporal_integrator>=0)
-    // 0=explicit trapezoidal predictor/corrector
-    // 1=Crank-Nicolson semi-implicit
-    // 2=explicit midpoint
-    // 3=multinomial diffusion
-    // 4=forward Euler  
     extern AMREX_GPU_MANAGED int reactDiff_diffusion_type;
-
-    // only used for split schemes (temporal_integrator>=0)
-    // 0=first-order (deterministic, tau leaping, CLE, or SSA)
-    // 1=second-order (determinisitc, tau leaping, or CLE only)
     extern AMREX_GPU_MANAGED int reactDiff_reaction_type;
-
-    // only used for midpoint diffusion schemes (split as well as unsplit)
-    // corrector formulation of noise
-    // 1 = K(nold) * W1 + K(nold)         * W2
-    // 2 = K(nold) * W1 + K(npred)        * W2
-    // 3 = K(nold) * W1 + K(2*npred-nold) * W2
     extern AMREX_GPU_MANAGED int midpoint_stoch_flux_type;
-
-    // how to compute n on faces for stochastic weighting
-    // 1=arithmetic (with C0-Heaviside), 2=geometric, 3=harmonic
-    // 10=arithmetic average with discontinuous Heaviside function
-    // 11=arithmetic average with C1-smoothed Heaviside function
-    // 12=arithmetic average with C2-smoothed Heaviside function
     extern AMREX_GPU_MANAGED int avg_type;
-
-    // use the Einkemmer boundary condition fix (split schemes only)
     extern AMREX_GPU_MANAGED int inhomogeneous_bc_fix;
-
-    // volume multiplier (dv = product(dx(1:MAX_SPACEDIM))*volume_factor)
-    // only really intended for 3D since in 2D one can control the cell depth
     extern AMREX_GPU_MANAGED amrex::Real volume_factor;
-
-    // initial values to be used in init_n.f90
     extern AMREX_GPU_MANAGED Array2D<amrex::Real, 0, 2 ,0, MAX_SPECIES> n_init_in;
-
-    // initialize from model file
     extern AMREX_GPU_MANAGED int model_file_init;
-
-    // initialize with all number of molecules strictly integer
     extern AMREX_GPU_MANAGED int integer_populations;
-
-    // Fickian diffusion coeffs
     extern AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real, MAX_SPECIES> D_Fick;
-
-    // diffusion boundary stencil order
     extern AMREX_GPU_MANAGED int diffusion_stencil_order;
-
-    // implicit diffusion solve verbosity
     extern AMREX_GPU_MANAGED int diffusion_verbose;
-
-    // implicit diffusion solve bottom solver verbosity
     extern AMREX_GPU_MANAGED int diffusion_bottom_verbose;
-
-    // relative eps for implicit diffusion solve
     extern AMREX_GPU_MANAGED amrex::Real implicit_diffusion_rel_eps;
-
-    // absolute eps for implicit diffusion solve
     extern AMREX_GPU_MANAGED amrex::Real implicit_diffusion_abs_eps;
-
 }