From 1dcdf2805a6c2a22975d51ac4d12fa3481613bd3 Mon Sep 17 00:00:00 2001
From: prkkumar <prabhatams@gmail.com>
Date: Wed, 28 Aug 2024 16:12:37 -0700
Subject: [PATCH] input file and some minor tweaks

---
 .../Exec/inputs_sundials_mri                  | 28 +++++++++++++------
 .../Reaction-Diffusion/Source/main.cpp        | 28 ++++++++++++++++---
 2 files changed, 43 insertions(+), 13 deletions(-)

diff --git a/ExampleCodes/SUNDIALS/Reaction-Diffusion/Exec/inputs_sundials_mri b/ExampleCodes/SUNDIALS/Reaction-Diffusion/Exec/inputs_sundials_mri
index 17b382b5..96860e0d 100644
--- a/ExampleCodes/SUNDIALS/Reaction-Diffusion/Exec/inputs_sundials_mri
+++ b/ExampleCodes/SUNDIALS/Reaction-Diffusion/Exec/inputs_sundials_mri
@@ -2,14 +2,18 @@
 n_cell = 32
 max_grid_size = 16
 
-nsteps = 1000
-plot_int = 100
+nsteps = 12
+plot_int = 2
 
-dt = 1.e-5
+dt = 8.5e-3
 
 # To replicate heat equation
 diffusion_coef = 1.0
-reaction_coef = 1.0
+reaction_coef = 1.e-4
+
+# MRI parameters
+use_MRI = true
+fast_dt_ratio = 0.02
 
 # Use adaptive time stepping and set integrator relative and absolute tolerances
 # adapt_dt = true
@@ -36,11 +40,17 @@ integration.type = SUNDIALS
 #
 # Optionally select a specific SUNDIALS method by name, see the SUNDIALS
 # documentation for the supported method names
+integration.sundials.type = EX-MRI
+integration.sundials.fast_type = ERK
 
-# Use the SUNDIALS default method for the chosen type (fixed or adaptive step sizes)
-integration.sundials.type = EX-MRI # all the types listed above work with respective default methods
+## *** Select a specific SUNDIALS ERK method ***
+#integration.sundials.method = ARKODE_BOGACKI_SHAMPINE_4_2_3
+#
+## *** Select a specific SUNDIALS ImEx method ***
+#integration.sundials.method_i = ARKODE_ARK2_DIRK_3_1_2
+#integration.sundials.method_e = ARKODE_ARK2_ERK_3_1_2
 
-# The following combination of type and method works as well
-# integration.sundials.type = EX-MRI
-# integration.sundials.method = ARKODE_MRI_GARK_FORWARD_EULER
+# *** Select a specific SUNDIALS MRI method ***
+integration.sundials.method = ARKODE_MIS_KW3
+integration.sundials.fast_method = ARKODE_KNOTH_WOLKE_3_3
 
diff --git a/ExampleCodes/SUNDIALS/Reaction-Diffusion/Source/main.cpp b/ExampleCodes/SUNDIALS/Reaction-Diffusion/Source/main.cpp
index 3a77eee4..f4f6552b 100644
--- a/ExampleCodes/SUNDIALS/Reaction-Diffusion/Source/main.cpp
+++ b/ExampleCodes/SUNDIALS/Reaction-Diffusion/Source/main.cpp
@@ -40,6 +40,10 @@ void main_main ()
     // use adaptive time step (dt used to set output times)
     bool adapt_dt = false;
 
+    // use MRI
+    bool use_MRI = false;
+    Real fast_dt_ratio = 0.1;
+
     // adaptive time step relative and absolute tolerances
     Real reltol = 1.0e-4;
     Real abstol = 1.0e-9;
@@ -73,6 +77,10 @@ void main_main ()
         // use adaptive step sizes
         pp.query("adapt_dt",adapt_dt);
 
+        // use MRI
+        pp.query("use_MRI",use_MRI);
+        pp.query("fast_dt_ratio",fast_dt_ratio);
+
         // adaptive step tolerances
         pp.query("reltol",reltol);
         pp.query("abstol",abstol);
@@ -177,7 +185,6 @@ void main_main ()
         Real reaction_coef = 0.0;
         Real diffusion_coef = 1.0;
         ParmParse pp;
-        pp.query("reaction_coef",reaction_coef);
         pp.query("diffusion_coef",diffusion_coef);
 
         for ( MFIter mfi(phi_data); mfi.isValid(); ++mfi )
@@ -226,7 +233,16 @@ void main_main ()
             // fill the right-hand-side for phi
             amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
             {
-                phi_rhs_array(i,j,k) = -1.*reaction_coef*phi_array(i,j,k);
+                if (use_MRI) {
+                    phi_rhs_array(i,j,k) = -1.*reaction_coef*phi_array(i,j,k);
+                } else {
+                    phi_rhs_array(i,j,k) = diffusion_coef*( (phi_array(i+1,j,k) - 2.*phi_array(i,j,k) + phi_array(i-1,j,k)) / (dx[0]*dx[0])
+                                                           +(phi_array(i,j+1,k) - 2.*phi_array(i,j,k) + phi_array(i,j-1,k)) / (dx[1]*dx[1])
+#if (AMREX_SPACEDIM == 3)
+                                                           +(phi_array(i,j,k+1) - 2.*phi_array(i,j,k) + phi_array(i,j,k-1)) / (dx[2]*dx[2]) )
+#endif
+                                         -1.*reaction_coef*phi_array(i,j,k);
+                }           
             });
         }
     };
@@ -234,7 +250,9 @@ void main_main ()
 
     TimeIntegrator<MultiFab> integrator(phi, time);
     integrator.set_rhs(rhs_function);
-    integrator.set_fast_rhs(rhs_fast_function);
+    if (use_MRI) {
+        integrator.set_fast_rhs(rhs_fast_function);
+    }
     if (adapt_dt) {
         integrator.set_adaptive_step();
         integrator.set_tolerances(reltol, abstol);
@@ -242,7 +260,9 @@ void main_main ()
         integrator.set_time_step(dt);
     }
 
-    integrator.set_fast_time_step(0.1*dt);
+    if (use_MRI) {
+        integrator.set_fast_time_step(fast_dt_ratio*dt);
+    }
 
     for (int step = 1; step <= nsteps; ++step)
     {