diff --git a/src/driver/driver.cpp b/src/driver/driver.cpp index 3ae676de..4b059c4f 100644 --- a/src/driver/driver.cpp +++ b/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__ diff --git a/src/driver/driver.hpp b/src/driver/driver.hpp index bd0ce35f..7fcc7448 100644 --- a/src/driver/driver.hpp +++ b/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; diff --git a/src/ion-neutral/ion-neutral.hpp b/src/ion-neutral/ion-neutral.hpp index 92c9af6d..b09ad325 100644 --- a/src/ion-neutral/ion-neutral.hpp +++ b/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 { diff --git a/src/ion-neutral/ion-neutral_tasks.cpp b/src/ion-neutral/ion-neutral_tasks.cpp index 2ae498b1..5f682cf8 100644 --- a/src/ion-neutral/ion-neutral_tasks.cpp +++ b/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),