fast rhs fun

AMReX-Codes · Jul 17, 2024 · ffa94a0 · ffa94a0
1 parent 21c5d47
commit ffa94a0
Showing 1 changed file with 38 additions and 0 deletions.
diff --git a/ExampleCodes/SUNDIALS/Reaction-Diffusion/Source/main.cpp b/ExampleCodes/SUNDIALS/Reaction-Diffusion/Source/main.cpp
@@ -202,16 +202,54 @@ void main_main ()
         }
     };
 
+    auto rhs_fast_function = [&](MultiFab& S_rhs, const MultiFab& S_data,
+                            const Real /* time */) {
+
+        // loop over boxes
+        auto& phi_data = S_data;
+        auto& phi_rhs  = S_rhs;
+
+        //Reaction and Diffusion Coefficients
+        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 )
+        {
+            const Box& bx = mfi.validbox();
+
+            const Array4<const Real>& phi_array = phi_data.array(mfi);
+            const Array4<Real>& phi_rhs_array = phi_rhs.array(mfi);
+
+            // 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) = 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
+            });
+        }
+    };
+
     TimeIntegrator<MultiFab> integrator(phi, time);
     integrator.set_pre_rhs_action(pre_rhs_function);
     integrator.set_rhs(rhs_function);
+    integrator.set_fast_rhs(rhs_fast_function);
     if (adapt_dt) {
         integrator.set_adaptive_step();
         integrator.set_tolerances(reltol, abstol);
     } else {
         integrator.set_time_step(dt);
     }
 
+    // This sets the ratio of slow timestep size to fast timestep size as an integer,
+    // or equivalently, the number of fast timesteps per slow timestep.
+    integrator.set_slow_fast_timestep_ratio(10);
+
     for (int step = 1; step <= nsteps; ++step)
     {
         // Set time to evolve to