setup preliminaries

exec/reactDiff/inputs_2d

@@ -0,0 +1,28 @@
+  # Problem specification
+  prob_lo = 0.0 0.0       # physical lo coordinate
+  prob_hi = 1.0 1.0       # physical hi coordinate
+
+  # number of cells in domain
+  n_cells = 32 32
+  # max number of cells in a box
+  max_grid_size = 16 16
+
+  # Time-step control
+  fixed_dt = 0.1
+
+  # Controls for number of steps between actions
+  max_step = 10
+  plot_int = 1
+
+  seed = 1
+
+  nspecies = 2
+
+  # Boundary conditions
+  # ----------------------
+  # BC specifications:
+  # -1 = periodic
+  bc_spec_lo = -1 -1
+  bc_spec_hi = -1 -1
+
+  stats_int = -1

exec/reactDiff/inputs_3d

@@ -0,0 +1,28 @@
+  # Problem specification
+  prob_lo = 0.0 0.0 0.0       # physical lo coordinate
+  prob_hi = 1.0 1.0 1.0       # physical hi coordinate
+
+  # number of cells in domain
+  n_cells = 32 32 32
+  # max number of cells in a box
+  max_grid_size = 16 16 16
+
+  # Time-step control
+  fixed_dt = 0.1
+
+  # Controls for number of steps between actions
+  max_step = 10
+  plot_int = 1
+
+  seed = 1
+
+  nspecies = 2
+
+  # Boundary conditions
+  # ----------------------
+  # BC specifications:
+  # -1 = periodic
+  bc_spec_lo = -1 -1 -1
+  bc_spec_hi = -1 -1 -1
+
+  stats_int = -1

exec/reactDiff/main_driver.cpp

@@ -1,5 +1,7 @@
 
 #include "common_functions.H"
+#include "chemistry_functions.H"
+#include "reactDiff_functions.H"
 
 #include <AMReX_VisMF.H>
 #include <AMReX_PlotFileUtil.H>
@@ -10,6 +12,8 @@
 
 using namespace std::chrono;
 using namespace amrex;
+using namespace common;
+using namespace chemistry;
 
 // argv contains the name of the inputs file entered at the command line
 void main_driver(const char* argv)
@@ -24,6 +28,8 @@ void main_driver(const char* argv)
 
     // Initialize variables in namespaces
     InitializeCommonNamespace();
+    InitializeChemistryNamespace();
+    InitializeReactDiffNamespace();
 
     // is the problem periodic?
     Vector<int> is_periodic(AMREX_SPACEDIM,0);  // set to 0 (not periodic) by default
@@ -42,14 +48,7 @@ void main_driver(const char* argv)
     Box domain(dom_lo, dom_hi);
 
     Geometry geom(domain,&real_box,CoordSys::cartesian,is_periodic.data());
-
-    // BoxArray
-    BoxArray ba;
-
-    // how boxes are distrubuted among MPI processes
-    DistributionMapping dmap;
 
-    Real dt = fixed_dt;
     const Real* dx = geom.CellSize();
 
     /////////////////////////////////////////
@@ -76,25 +75,90 @@ void main_driver(const char* argv)
         }
 
     }
+
+    BoxArray ba;
+    DistributionMapping dmap;
 
     int step_start;
     amrex::Real time;
 
-    // Initialize the boxarray "ba" from the single box "bx"
-    ba.define(domain);
+    if (restart < 0) {
+
+        step_start = 1;
+        time = 0.;
+
+        // Initialize the boxarray "ba" from the single box "bx"
+        ba.define(domain);
+
+        // Break up boxarray "ba" into chunks no larger than "max_grid_size" along a direction
+        // note we are converting "Vector<int> max_grid_size" to an IntVect
+        ba.maxSize(IntVect(max_grid_size));
 
-    // Break up boxarray "ba" into chunks no larger than "max_grid_size" along a direction
-    // note we are converting "Vector<int> max_grid_size" to an IntVect
-    ba.maxSize(IntVect(max_grid_size));
+        dmap.define(ba);
+
+    } else {
+
+        // checkpoint restart
+
+    }
 
-    dmap.define(ba);
+    MultiFab n_old(ba,dmap,nspecies,1);
+    MultiFab n_new(ba,dmap,nspecies,1);
+
+    if (model_file_init) {
+        Abort("model_file_init not supported yet");
+    } else {
+        // Initialize n
+        // Init();
+        n_old.setVal(0.);
+    }
 
+    if (std::abs(initial_variance_mass) > 0.) {
+        Abort("initial_variance_mass not supported yet");
+        // add_init_n_fluctuations()
+    }
 
+    Real dt;
+    if (fixed_dt > 0.) {
+        dt = fixed_dt;
+        Print() << "Setting dt using fixed_dt = " << dt << std::endl;
+    } else {
+        Real D_Fick_max = 0.;
+        for (int i=0; i<nspecies; ++i ) {
+            D_Fick_max = std::max(D_Fick_max,D_Fick[i]);
+        }
+        Real dx_min = dx[0];
+        for (int i=1; i<AMREX_SPACEDIM; ++i) {
+            dx_min = std::min(dx_min,dx[i]);
+        }
+        dt = cfl * dx_min / (2. * AMREX_SPACEDIM * D_Fick_max);
+        Print() << "Setting dt using explicit diffusion cfl condition = " << dt << std::endl;
+    }
 
+    if (inhomogeneous_bc_fix == 1 && temporal_integrator > 0) {
+        Abort("comput_n_steady not supported for inhomogeneous_bc_fix == 1 && temporal_integrator > 0 yet");
+        // compute_n_steady()
+    }
 
+    if (temporal_integrator < 0) { // unsplit schemes
+        // Donev: The code will work for a single cell also but may not be the most efficient, so issue warning:
+        if (n_cells[0] == 1 && n_cells[1] == 1) {
+            Print() << "WARNING in advance_reaction_diffusion: use splitting based schemes (temporal_integrator>=0) for single cell" << std::endl;
+        }
+        if (nreaction < 1) {
+            Print() << "WARNING in advance_reaction_diffusion: use splitting based schemes (temporal_integrator>=0) for diffusion only" << std::endl;
+        }
+    }
 
+    if (stats_int > 0) {
+        Abort("Structure factor not implemented yet");
+    }
+
+    int istep = (restart < 0) ? 0 : restart;
+    WritePlotFile(istep,time,geom,n_old);
 
 
+
     ///////////////////////////////////////////
 
     // time step loop

src_reactDiff/Make.package

@@ -1,3 +1,4 @@
 CEXE_sources   += reactDiff_functions.cpp
+CEXE_sources   += WritePlotFile.cpp
 CEXE_headers   += reactDiff_functions.H
 CEXE_headers   += reactDiff_namespace.H

src_reactDiff/WritePlotFile.cpp

@@ -0,0 +1,32 @@
+#include "AMReX_PlotFileUtil.H"
+#include "reactDiff_functions.H"
+
+void WritePlotFile(int step,
+                   const amrex::Real time,
+                   const amrex::Geometry& geom,
+		   const MultiFab& n_in)
+{
+
+    BL_PROFILE_VAR("WritePlotFile()",WritePlotFile);
+
+    std::string plotfilename = Concatenate(plot_base_name,step,7);
+
+    amrex::Print() << "Writing plotfile " << plotfilename << "\n";
+
+    BoxArray ba = n_in.boxArray();
+    DistributionMapping dmap = n_in.DistributionMap();
+
+    Vector<std::string> varNames(nspecies);
+
+    // keep a counter for plotfile variables
+    int cnt = 0;
+
+    for (int i=0; i<nspecies; ++i) {
+        std::string x = "n";
+        x += char(49+i);
+        varNames[cnt++] = x;
+    }
+
+    WriteSingleLevelPlotfile(plotfilename,n_in,varNames,geom,time,step);
+
+}

src_reactDiff/reactDiff_functions.H

@@ -13,4 +13,8 @@ using namespace common;
 
 void InitializeReactDiffNamespace();
 
+void WritePlotFile(int step,
+                   const amrex::Real time,
+                   const amrex::Geometry& geom,
+		   const MultiFab& n_in);
 #endif

src_reactDiff/reactDiff_functions.cpp

@@ -47,11 +47,14 @@ AMREX_GPU_MANAGED amrex::Real reactDiff::volume_factor;
 // initial values to be used in init_n.f90
 AMREX_GPU_MANAGED Array2D<amrex::Real, 0, 2 ,0, MAX_SPECIES> reactDiff::n_init_in;
 
+// initialize from model file
+AMREX_GPU_MANAGED int reactDiff::model_file_init;
+
 // initialize with all number of molecules strictly integer
 AMREX_GPU_MANAGED int reactDiff::integer_populations;
 
 // Fickian diffusion coeffs
-AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real, MAX_SPECIES> reactDiff::D_fick;
+AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real, MAX_SPECIES> reactDiff::D_Fick;
 
 // diffusion boundary stencil order
 AMREX_GPU_MANAGED int reactDiff::diffusion_stencil_order;
@@ -110,12 +113,15 @@ void InitializeReactDiffNamespace()
         }
     }
 
+    model_file_init = 0;
+    pp.query("model_file_init",model_file_init);
+
     integer_populations = 0;
     pp.query("integer_populations",integer_populations);
 
-    if (pp.queryarr("D_fick",temp)) {
+    if (pp.queryarr("D_Fick",temp)) {
         for (int i=0; i<nspecies; ++i) {
-            D_fick[i] = temp[i];
+            D_Fick[i] = temp[i];
         }
     }
 

src_reactDiff/reactDiff_namespace.H

@@ -46,11 +46,14 @@ namespace reactDiff {
     // 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;
+    extern AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real, MAX_SPECIES> D_Fick;
 
     // diffusion boundary stencil order
     extern AMREX_GPU_MANAGED int diffusion_stencil_order;