Implicit chemisty implementation

Source/mflo.H

@@ -76,6 +76,14 @@ public:
     ErrorEst(int lev, amrex::TagBoxArray& tags, amrex::Real time, int ngrow)
         override;
 
+    // advance a single level for a single time step, updates flux registers
+    void Advance_coupled_strang(
+        int lev,
+        amrex::Real time,
+        amrex::Real dt_lev,
+        int iteration,
+        int ncycle,bool only_flow=false);
+
     // advance a single level for a single time step, updates flux registers
     void Advance_coupled(
         int lev,
@@ -87,6 +95,9 @@ public:
     void Advance_chemistry(int lev, 
             Real time, Real dt_lev);
 
+    void Advance_chemistry_implicit(int lev,
+            Real time, Real dt_lev);
+
     // compute dt from CFL considerations
     Real EstTimeStep(int lev);
 

Source/mflo.cpp

@@ -293,6 +293,16 @@ void mflo::ReadParameters()
         pp.query("dissfactor",dissfactor);
         pp.query("species_in_solid",species_in_solid);
     }
+    {
+        // Add default time integration options for the chemisty solver
+        // Note that this is explicit RK3, which does not require sundials
+        ParmParse pp("integration");
+        if( not pp.contains("type") ) {
+            pp.add("type", "RungeKutta");
+            pp.add("rk.type", 3);
+            Print() << "Defaulting to explicit RK3 scheme for chemistry solve" << std::endl;
+        }
+    }
 }
 
 // utility to copy in data from phi_old and/or phi_new into another multifab

Source/mflo_evolve.cpp

@@ -5,6 +5,7 @@
 #include <AMReX_PlotFileUtil.H>
 #include <AMReX_VisMF.H>
 #include <AMReX_PhysBCFunct.H>
+#include <AMReX_TimeIntegrator.H>
 #include <kernels_3d.H>
 
 #ifdef AMREX_MEM_PROFILING
@@ -355,7 +356,7 @@ void mflo::timeStep(int lev, Real time, int iteration, bool only_flow)
 
 
     // advance a single level for a single time step, updates flux registers
-    Advance_coupled(lev, time, dt[lev], iteration, nsubsteps[lev],only_flow);
+    Advance_coupled_strang(lev, time, dt[lev], iteration, nsubsteps[lev],only_flow);
 
     ++istep[lev];
 
@@ -385,6 +386,75 @@ void mflo::timeStep(int lev, Real time, int iteration, bool only_flow)
     }
 }
 
+void mflo::Advance_coupled_strang(int lev, Real time, Real dt_lev, int iteration, int ncycle,bool only_flow)
+{
+    constexpr int num_grow = 3;
+    std::swap(phi_old[lev], phi_new[lev]); // old becomes new and new becomes
+                                           // old
+    t_old[lev] = t_new[lev];        // old time is now current time (time)
+    t_new[lev] += dt_lev;           // new time is ahead
+    MultiFab& S_new = phi_new[lev]; // this is the old value, beware!
+    MultiFab& S_old = phi_old[lev]; // current value
+
+    //RK3 TVD scheme
+    // State with ghost cells
+    MultiFab Sborder(grids[lev], dmap[lev], S_new.nComp(), num_grow);
+    // source term
+    MultiFab dsdt_flow(grids[lev], dmap[lev], S_new.nComp(), 0);
+
+    if (do_reflux)
+    {
+        if (flux_reg[lev + 1])
+        {
+            flux_reg[lev+1]->setVal(Real(0.0));
+        }
+    }
+
+
+    // stage 1
+    // time is current time which is t_old
+    FillPatch(lev, time, Sborder, 0, Sborder.nComp());
+    // compute dsdt for 1/2 timestep
+    update_cutcell_data(lev, num_grow, Sborder, time, dt_lev);
+    compute_dsdt_flow(lev, num_grow, Sborder, dsdt_flow, time, dt_lev, sixth, true);
+    // S_new=S_old+dt*dsdt //sold is the current value
+    MultiFab::LinComb(S_new, one, Sborder, 0, dt_lev, dsdt_flow, 0, 0, S_new.nComp(), 0);
+
+
+    // stage 2
+    // time+dt_lev lets me pick S_new for sborder
+    FillPatch(lev, time + dt_lev, Sborder, 0, Sborder.nComp());
+    update_cutcell_data(lev,num_grow,Sborder,time,dt_lev);
+    compute_dsdt_flow(lev, num_grow, Sborder, dsdt_flow, time + dt_lev, dt_lev, sixth,true);
+
+    // S_new=3/4 S_old+1/4 S_new + 1/4 dt*dsdt
+    // S_new = S_new + dt*dsdt
+    // S_new = S_new + 3*S_old
+    // S_new =S_new/4
+    MultiFab::Saxpy(S_new, dt_lev, dsdt_flow, 0, 0, S_new.nComp(), 0);
+    MultiFab::Saxpy(S_new, Real(3.0), S_old, 0, 0, S_new.nComp(), 0);
+    S_new.mult(fourth);
+
+    // stage 3
+    // time+dt_lev lets me pick S_new for sborder
+    FillPatch(lev, time + dt_lev, Sborder, 0, Sborder.nComp());
+    update_cutcell_data(lev,num_grow,Sborder,time,dt_lev);
+    compute_dsdt_flow(lev, num_grow, Sborder, dsdt_flow, time + dt_lev, dt_lev, two3rd,true);
+    // S_new=1/3 S_old+2/3 S_new + 2/3 dt*dsdt
+    // S_new = S_new + dt*dsdt
+    // S_new = S_new + 2*S_old
+    // S_new =S_new/3
+    MultiFab::Saxpy(S_new, dt_lev, dsdt_flow, 0, 0, S_new.nComp(), 0);
+    MultiFab::Xpay(S_new, two, S_old, 0, 0, S_new.nComp(), 0);
+    S_new.mult(third);
+
+    //Advance chemistry here
+    if(!only_flow)
+    {
+        Advance_chemistry_implicit(lev, time, dt_lev);
+    }
+}
+
 void mflo::Advance_coupled(int lev, Real time, Real dt_lev, int iteration, int ncycle,bool only_flow)
 {
     constexpr int num_grow = 3;
@@ -455,6 +525,8 @@ void mflo::Advance_coupled(int lev, Real time, Real dt_lev, int iteration, int n
     // time+dt_lev lets me pick S_new for sborder
     FillPatch(lev, time + dt_lev, Sborder, 0, Sborder.nComp());
     update_cutcell_data(lev,num_grow,Sborder,time,dt_lev);
+    //NOTE: Using time + dt_lev might be incorrect here, need to verify using beta from the RK3 coefficients
+    //According to one source, this might need to be t + 0.5*dt_lev
     compute_dsdt_flow(lev, num_grow, Sborder, dsdt_flow, time + dt_lev, dt_lev, two3rd,true);
     if(!only_flow)
     {
@@ -545,6 +617,45 @@ void mflo::Advance_chemistry(int lev, Real time, Real dt_lev)
 
 }
 
+/*
+Advances the chemisty state from time -> time + dt_lev
+This is done using TimeIntegrator from AMReX, which can
+be either implicit or explicit depending on the parameters
+in the input file. Implicit currently requires an existing
+SUNDIALS installation, with USE_SUNDIALS=TRUE at compile time,
+along with setting the following in the input file
+
+Explicit (no sundials):
+   integration.type = RungeKutta
+   integration.rk.type = 3 (for 3rd order)
+Implicit (with sundials):
+   integration.type = SUNDIALS
+   integration.sundials.strategy = CVODE
+ */
+void mflo::Advance_chemistry_implicit(int lev, Real time, Real dt_lev)
+{
+    constexpr int num_grow = 3;
+    std::swap(phi_old[lev], phi_new[lev]); // old becomes new and new becomes old
+    MultiFab& S_new = phi_new[lev]; // old value
+    MultiFab& S_old = phi_old[lev]; // current value
+
+    auto rhs_function = [&] ( Vector<MultiFab> & dSdt_vec, const Vector<MultiFab>& S_vec, const Real time) {
+        auto & dSdt = dSdt_vec[0];
+        MultiFab S(S_vec[0], amrex::make_alias, 0, S_vec[0].nComp());
+        compute_dsdt_chemistry(lev, num_grow, S, dSdt, time);
+    };
+    Vector<MultiFab> state_old, state_new;
+    // This term has the current state
+    state_old.push_back(MultiFab(S_old, amrex::make_alias, 0, S_new.nComp()));
+    // This is where the integrator puts the new state, hence aliased to S_new
+    state_new.push_back(MultiFab(S_new, amrex::make_alias, 0, S_new.nComp()));
+    // Define the integrator
+    TimeIntegrator<Vector<MultiFab>> integrator(state_old);
+    integrator.set_rhs(rhs_function);
+    // Advance from time to time + dt_lev
+    integrator.advance(state_old, state_new, time, dt_lev); //S_new/phi_new should have the new state
+}
+
 void mflo::compute_dsdt_chemistry(
     int lev,
     const int num_grow,

models/blockCatalyst1d/GNUmakefile

@@ -4,15 +4,34 @@ PROFILE    = FALSE
 DEBUG      = TRUE
 DEBUG      = FALSE
 
-DIM       = 3
+DIM        = 3
 
-COMP	   = gnu
+COMP       = gnu
 
 USE_MPI    = TRUE
-USE_OMP    = FALSE 
+USE_OMP    = FALSE
 USE_CUDA   = FALSE
 
-Bpack   := ./Make.package 
-Blocs   := . 
+Bpack   := ./Make.package
+Blocs   := .
 
 include ../Make.mflo
+
+ifeq ($(USE_SUNDIALS),TRUE)
+# NOTE: SUNDIALS_ROOT must point to the directory where sundials is installed
+# A good check is to see if $(SUNDIALS_ROOT)/lib has a bunch of libsundials_ files
+# To run with sundials, enabled, please compile with USE_SUNDIALS = TRUE
+SUNDIALS_ROOT ?= $(TOP)/../../../../../sundials/instdir
+SUNDIALS_LIB_DIR ?= $(SUNDIALS_ROOT)/lib
+
+USE_CVODE_LIBS ?= TRUE
+USE_ARKODE_LIBS ?= TRUE
+
+DEFINES += -DAMREX_USE_SUNDIALS
+INCLUDE_LOCATIONS += $(SUNDIALS_ROOT)/include
+LIBRARY_LOCATIONS += $(SUNDIALS_LIB_DIR)
+
+LIBRARIES += -L$(SUNDIALS_LIB_DIR) -lsundials_cvode
+LIBRARIES += -L$(SUNDIALS_LIB_DIR) -lsundials_arkode
+LIBRARIES += -L$(SUNDIALS_LIB_DIR) -lsundials_nvecmanyvector
+endif

models/blockCatalyst1d/inputs_x

@@ -52,3 +52,10 @@ amr.plot_int   = 10000     # number of timesteps between plot files
 # CHECKPOINT
 amr.chk_file = chk      # root name of checkpoint file
 amr.chk_int  = -1       # number of timesteps between checkpoint files
+
+# INTEGRATION options (for chemistry)
+## Uncomment below two lines to enable implicit chemistry
+# integration.type = SUNDIALS
+# integration.sundials.strategy = CVODE
+## Note that this requires a working sundials installation and
+## the code must be compiled with USE_SUNDIALS=TRUE, see GNUmakefile for additional info

models/blockCatalyst1d/inputs_y

@@ -52,3 +52,10 @@ amr.plot_int   = 10000     # number of timesteps between plot files
 # CHECKPOINT
 amr.chk_file = chk      # root name of checkpoint file
 amr.chk_int  = -1       # number of timesteps between checkpoint files
+
+# INTEGRATION options (for chemistry)
+## Uncomment below two lines to enable implicit chemistry
+# integration.type = SUNDIALS
+# integration.sundials.strategy = CVODE
+## Note that this requires a working sundials installation and
+## the code must be compiled with USE_SUNDIALS=TRUE, see GNUmakefile for additional info

models/blockCatalyst1d/inputs_z

@@ -52,3 +52,10 @@ amr.plot_int   = 10000     # number of timesteps between plot files
 # CHECKPOINT
 amr.chk_file = chk      # root name of checkpoint file
 amr.chk_int  = -1       # number of timesteps between checkpoint files
+
+# INTEGRATION options (for chemistry)
+## Uncomment below two lines to enable implicit chemistry
+# integration.type = SUNDIALS
+# integration.sundials.strategy = CVODE
+## Note that this requires a working sundials installation and
+## the code must be compiled with USE_SUNDIALS=TRUE, see GNUmakefile for additional info

models/blockCatalyst1d/runcase.sh

@@ -1,13 +1,20 @@
 #!/bin/bash
+# USAGE example: sh runcase.sh mpiexec 8
+# run in x direction
 $1 -n $2 ./*.ex inputs_x
 finfile=$(ls -d plt?????/ | tail -n 1)
-mv ${finfile} finalplt_x
-. ../clean.sh
+rm -rf finalplt_x #Remove existing finalplt_x directory
+mv ${finfile} finalplt_x #Move the last plt file for plotting
+. ./clean.sh #Remove plt/chk files
+
+# Repeat for other two directions
 $1 -n $2 ./*.ex inputs_y
 finfile=$(ls -d plt?????/ | tail -n 1)
+rm -rf finalplt_y
 mv ${finfile} finalplt_y
-. ../clean.sh
+. ./clean.sh
 $1 -n $2 ./*.ex inputs_z
 finfile=$(ls -d plt?????/ | tail -n 1)
+rm -rf finalplt_z
 mv ${finfile} finalplt_z
-. ../clean.sh
+. ./clean.sh