New IMEX integrator

src/driver/driver.cpp

@@ -165,8 +165,8 @@ Driver::Driver(ParameterInput *pin, Mesh *pmesh, Real wtlim, Kokkos::Timer* ptim
       nimp_stages = 3;
       nexp_stages = 2;
       cfl_limit = 1.0;
-      gam0[0] = 0.0;
-      gam1[0] = 1.0;
+      gam0[0] = 1.0;
+      gam1[0] = 0.0;
       beta[0] = 1.0;
 
       gam0[1] = 0.5;
@@ -185,6 +185,48 @@ Driver::Driver(ParameterInput *pin, Mesh *pmesh, Real wtlim, Kokkos::Timer* ptim
       a_twid[2][1] = 0.25;
       a_twid[2][2] = 0.25;
       a_impl = 0.5;
+    } else if (integrator == "imex2+") {
+      // IMEX(4,3,2): Krapp et al. (2024, arXiv:2310.04435), Eq.30.
+      // three-stage explicit, four-stage implicit, second-order ImEx
+      // two implicit stages added, adapting Athenak's overall architecture
+      // Note explicit steps may not reduce to RK2 based on the parameters chosen
+      nimp_stages = 4;
+      nexp_stages = 3;
+      cfl_limit = 1.0;
+      gamma = 1.707106781186547;   //1+1/sqrt(2)
+      gam0[0] = 1.0;
+      gam1[0] = 0.0;
+      beta[0] = gamma;
+
+      gam0[1] = (2.0*gamma-1.0)/(2.0*gamma*gamma);
+      gam1[1] = (1.0-(2.0*gamma-1.0)/(2.0*gamma*gamma));
+      beta[1] = 1.0/(2.0*gamma);
+
+      gam0[2] = 1.0;
+      gam1[2] = 0.0;
+      beta[2] = 0.0;
+
+      a_twid[0][0] = 0.0;
+      a_twid[0][1] = 0.0;
+      a_twid[0][2] = 0.0;
+      a_twid[0][3] = 0.0;
+
+      a_twid[1][0] = 0.0;
+      a_twid[1][1] = 0.0;
+      a_twid[1][2] = 0.0;
+      a_twid[1][3] = 0.0;
+
+      a_twid[2][0] = 0.0;
+      a_twid[2][1] = 0.0;
+      a_twid[2][2] = (1.0-2.0*gamma*gamma)/2.0/gamma;
+      a_twid[2][3] = 0.0;
+
+      a_twid[3][0] = 0.0;
+      a_twid[3][1] = 0.0;
+      a_twid[3][2] = 0.0;
+      a_twid[3][3] = 0.0;
+
+      a_impl = gamma;
     } else if (integrator == "imex3") {
       // IMEX-SSP3(4,3,3): Pareschi & Russo (2005) Table VI.
       // three-stage explicit, four-stage implicit, third-order ImEx
@@ -227,33 +269,6 @@ Driver::Driver(ParameterInput *pin, Mesh *pmesh, Real wtlim, Kokkos::Timer* ptim
       a_twid[3][2] = (4.0*(b + e + a) - 1.0)/6.0;
       a_twid[3][3] = 2.0*(1.0 - a)/3.0;
       a_impl = a;
-    } else if (integrator == "imex+") {
-      // IMEX(2,3,2): Krapp et al. (2024, arXiv:2310.04435), Eq.30.
-      // three-stage explicit, two-stage implicit, second-order ImEx
-      // Note explicit steps may not reduce to RK2 based on the parameters chosen
-      nimp_stages = 2;
-      nexp_stages = 3;
-      cfl_limit = 1.0;
-      Real gamma = 1.7071067811865475;   // 1 + 1/sqrt(2)
-      gam0[0] = 0.0;
-      gam1[0] = 1.0;
-      beta[0] = gamma;
-
-      gam0[1] = (2.0*gamma-1.0)/(2.0*gamma*gamma);
-      gam1[1] = 1.0-(2.0*gamma-1.0)/(2.0*gamma*gamma);
-      beta[1] = 1.0/(2.0*gamma);
-
-      gam0[2] = 1.0;
-      gam1[2] = 0.0;
-      beta[2] = 0.0;
-
-      a_twid[0][0] = (1.0-2*gamma*gamma)/(2.0*gamma);
-      a_twid[0][1] = 0.0;
-
-      a_twid[1][0] = 0.0;
-      a_twid[1][1] = 0.0;
-
-      a_impl = gamma;
     // Error, unrecognized integrator name.
     } else {
       std::cout << "### FATAL ERROR in " << __FILE__ << " at line " << __LINE__

src/driver/driver.hpp

@@ -44,6 +44,7 @@ class Driver {
   Real delta[4];                   // weights for updating the intermediate stage (u1)
   Real a_twid[4][4], a_impl;       // matrix elements for implicit stages in ImEx
   Real cfl_limit;                  // maximum CFL number for integrator
+  Real gamma;                      // gamma value for the IMEX_new integrator
   Kokkos::Timer* pwall_clock_;     // timer for tracking the wall clock
   Real wall_time;
 

src/ion-neutral/ion-neutral.hpp

@@ -57,6 +57,8 @@ struct IonNeutralTaskIDs {
   TaskID n_newdt;
   TaskID i_clear;
   TaskID n_clear;
+  TaskID i_srctrms;
+  TaskID n_srctrms;
 };
 
 namespace ion_neutral {

src/ion-neutral/ion-neutral_tasks.cpp

@@ -47,13 +47,15 @@ void IonNeutral::AssembleIonNeutralTasks(
   id.i_sendf  = tl["stagen"]->AddTask(&MHD::SendFlux, pmhd, id.i_flux);
   id.i_recvf  = tl["stagen"]->AddTask(&MHD::RecvFlux, pmhd, id.i_sendf);
   id.i_rkupdt = tl["stagen"]->AddTask(&MHD::RKUpdate, pmhd, id.i_recvf);
+  id.i_srctrms   = tl["stagen"]->AddTask(&MHD::MHDSrcTerms, pmhd, id.i_rkupdt);
 
-  id.n_flux   = tl["stagen"]->AddTask(&Hydro::Fluxes, phyd, id.i_rkupdt);
+  id.n_flux   = tl["stagen"]->AddTask(&Hydro::Fluxes, phyd, id.i_srctrms);
   id.n_sendf  = tl["stagen"]->AddTask(&Hydro::SendFlux, phyd, id.n_flux);
   id.n_recvf  = tl["stagen"]->AddTask(&Hydro::RecvFlux, phyd, id.n_sendf);
   id.n_rkupdt = tl["stagen"]->AddTask(&Hydro::RKUpdate, phyd, id.n_recvf);
+  id.n_srctrms   = tl["stagen"]->AddTask(&Hydro::HydroSrcTerms, phyd, id.n_rkupdt);
 
-  id.impl     = tl["stagen"]->AddTask(&IonNeutral::ImpRKUpdate, this, id.n_rkupdt);
+  id.impl     = tl["stagen"]->AddTask(&IonNeutral::ImpRKUpdate, this, id.n_srctrms);
   id.i_restu  = tl["stagen"]->AddTask(&MHD::RestrictU, pmhd, id.impl);
   id.n_restu  = tl["stagen"]->AddTask(&Hydro::RestrictU, phyd, id.i_restu);
 
@@ -185,9 +187,24 @@ TaskStatus IonNeutral::ImpRKUpdate(Driver *pdriver, int estage) {
   // equations for ion-neutral drag.
   // Only required for istage = (1,2,3,[4])
   if (estage < pdriver->nexp_stages) {
-    Real gamma_adt = drag_coeff*(pdriver->a_impl)*(pmy_pack->pmesh->dt);
-    Real xi_adt = ionization_coeff*(pdriver->a_impl)*(pmy_pack->pmesh->dt);
-    Real alpha_adt = recombination_coeff*(pdriver->a_impl)*(pmy_pack->pmesh->dt);
+    Real gamma_adt;
+    Real xi_adt;
+    Real alpha_adt;
+
+    // Condition to set gamma_adt, xi_adt, and alpha_adt to zero
+    if (istage < 3 && pdriver->integrator == "imex2+") {
+      gamma_adt = 0.0;
+      xi_adt = 0.0;
+      alpha_adt = 0.0;
+    } else {
+      gamma_adt = drag_coeff * (pdriver->a_impl) * (pmy_pack->pmesh->dt);
+      xi_adt = ionization_coeff * (pdriver->a_impl) * (pmy_pack->pmesh->dt);
+      alpha_adt = recombination_coeff * (pdriver->a_impl) * (pmy_pack->pmesh->dt);
+    }
+
+    //Real gamma_adt = drag_coeff*(pdriver->a_impl)*(pmy_pack->pmesh->dt);
+    //Real xi_adt = ionization_coeff*(pdriver->a_impl)*(pmy_pack->pmesh->dt);
+    //Real alpha_adt = recombination_coeff*(pdriver->a_impl)*(pmy_pack->pmesh->dt)
     auto ui = pmhd->u0;
     auto un = phyd->u0;
     par_for("imex_imp",DevExeSpace(),0,nmb1,0,(n3-1),0,(n2-1),0,(n1-1),