From 664bd230828ef32b57f4ceef5390a3162a01a143 Mon Sep 17 00:00:00 2001 From: Lin Shih Date: Wed, 1 Apr 2026 10:27:01 +0100 Subject: [PATCH 01/12] add vth diagnose --- src/neutral_mixed.cxx | 20 ++++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 8715e2b88..20ba35380 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -680,8 +680,9 @@ void NeutralMixed::finally(const Options& state) { ///////////////////////////////////////////////////// // Neutral density TRACE("Neutral density"); - ddt(Nn) = -FV::Div_par_mod(Nn, Vn, sound_speed, - pf_adv_par_ylow); // Parallel advection + ddtN_par_advection = -FV::Div_par_mod(Nn, Vn, sound_speed, + pf_adv_par_ylow); + ddt(Nn) = ddtN_par_advection; // Parallel advection // Perpendicular diffusion if (nonorthogonal_operators) { @@ -1165,6 +1166,21 @@ void NeutralMixed::outputVars(Options& state) { {"long_name", name + " pressure source"}, {"species", name}, {"source", "neutral_mixed"}}); + // Lin add + set_with_attrs(state[std::string("V") + name + std::string("th_pf")], Vnth_pf, + {{"time_dimension", "t"}, + {"units", "m / s"}, + {"conversion", Cs0}, + {"standard_name", "velocity"}, + {"long_name", name + " thermal velocity of neutrals for advection"}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("V") + name + std::string("th_hf")], Vnth_hf, + {{"time_dimension", "t"}, + {"units", "m / s"}, + {"conversion", Cs0}, + {"standard_name", "velocity"}, + {"long_name", name + " thermal velocity of neutrals for heat flux"}, + {"source", "neutral_mixed"}}); /////////////////////////////////////////////////// // Parallel flow diagnostics From 9669d560f2676d246ca0a8602bc350ce53a34177 Mon Sep 17 00:00:00 2001 From: Lin Shih Date: Wed, 1 Apr 2026 11:32:08 +0100 Subject: [PATCH 02/12] add physcis terms --- include/neutral_mixed.hxx | 19 +++++ src/neutral_mixed.cxx | 170 ++++++++------------------------------ 2 files changed, 52 insertions(+), 137 deletions(-) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index d78dc050a..2ca455452 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -45,6 +45,9 @@ private: std::vector collision_names; ///< Collisions used for collisionality std::string diffusion_collisions_mode; ///< Collision selection, either afn or multispecies + Field3D Vnth_pf = 0.0; + Field3D Vnth_hf = 0.0; + Field3D nu; ///< Collisionality to use for diffusion Field3D nu_pseudo_mfp; ///< Pseudo-collision frequency based on mean free path Field3D nu_total; ///< Total collision frequency used for diffusion, including @@ -105,6 +108,22 @@ private: bool output_ddt; ///< Save time derivatives? bool diagnose; ///< Save additional diagnostics? + //Physics terms + Field3D ddtN_par_advection; + Field3D ddtN_perp_diffusion; + + Field3D ddtPn_par_advection; + Field3D ddtPn_work_done; + Field3D ddtPn_perp_advection; + Field3D ddtPn_par_conduction;:w + Field3D ddtPn_perp_conduction; + + Field3D ddtNVn_par_advection; + Field3D ddtNVn_pressure_gradient; + Field3D ddtNVn_perp_advection; + Field3D par_viscosity_source; + Field3D perp_viscosity_source; + // Flow diagnostics Field3D pf_adv_perp_xlow, pf_adv_perp_ylow, pf_adv_par_ylow; Field3D mf_adv_perp_xlow, mf_adv_perp_ylow, mf_adv_par_ylow; diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 20ba35380..99301945f 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -512,9 +512,6 @@ void NeutralMixed::finally(const Options& state) { // New formulation: 1D particle flow in 3D maxwellian, 1D heat flow in 3D Maxwellian // See Stangeby - Field3D Vnth_pf = 0.0; - Field3D Vnth_hf = 0.0; - if (legacy_thermal_speed) { Vnth_pf = sqrt(Tnlim / AA); Vnth_hf = 3.0 / 2.0 * Vnth_pf; @@ -686,10 +683,11 @@ void NeutralMixed::finally(const Options& state) { // Perpendicular diffusion if (nonorthogonal_operators) { - ddt(Nn) += - Div_a_Grad_perp_nonorthog(DnnNn, logPnlim, pf_adv_perp_xlow, pf_adv_perp_ylow); + ddtN_perp_diffusion = Div_a_Grad_perp_nonorthog(DnnNn, logPnlim, pf_adv_perp_xlow, pf_adv_perp_ylow); + ddt(Nn) += ddtN_perp_diffusion; } else { - ddt(Nn) += Div_a_Grad_perp_flows(DnnNn, logPnlim, pf_adv_perp_xlow, pf_adv_perp_ylow); + ddtN_perp_diffusion = Div_a_Grad_perp_flows(DnnNn, logPnlim, pf_adv_perp_xlow, pf_adv_perp_ylow); + ddt(Nn) += ddtN_perp_diffusion; } Sn = density_source; // Save for possible output @@ -702,20 +700,21 @@ void NeutralMixed::finally(const Options& state) { // Neutral pressure TRACE("Neutral pressure"); - ddt(Pn) = -(5. / 3) + ddtPn_par_advection = -(5. / 3) * FV::Div_par_mod( // Parallel advection Pn, Vn, sound_speed, ef_adv_par_ylow) - + (2. / 3) * Vn * Grad_par(Pn); // Work done + ddtPn_work_done = (2. / 3) * Vn * Grad_par(Pn); // Work done + ddt(Pn) = ddtPn_par_advection + ddtPn_work_done; // Perpendicular advection of pressure if (nonorthogonal_operators) { - ddt(Pn) += - (5. / 3) + ddtPn_perp_advection = (5. / 3) * Div_a_Grad_perp_nonorthog(DnnPn, logPnlim, ef_adv_perp_xlow, ef_adv_perp_ylow); + ddt(Pn) += ddtPn_perp_advection; } else { - ddt(Pn) += - (5. / 3) + ddtPn_perp_advection = (5. / 3) * Div_a_Grad_perp_flows(DnnPn, logPnlim, ef_adv_perp_xlow, ef_adv_perp_ylow); + ddt(Pn) += ddtPn_perp_advection; } // The factor here is 5/2 as we're advecting internal energy and pressure. @@ -724,20 +723,23 @@ void NeutralMixed::finally(const Options& state) { ef_adv_perp_ylow *= 5. / 2; if (neutral_conduction) { - ddt(Pn) += (2. / 3) + ddtPn_par_conduction = (2. / 3) * Div_par_K_Grad_par_mod(kappa_n_par, Tn, // Parallel conduction ef_cond_par_ylow, false); // No conduction through target boundary + ddt(Pn) += ddtPn_par_conduction; // Perpendicular conduction if (nonorthogonal_operators) { - ddt(Pn) += (2. / 3) + ddtPn_perp_conduction = (2. / 3) * Div_a_Grad_perp_nonorthog(kappa_n_perp, Tn, ef_cond_perp_xlow, ef_cond_perp_ylow); + ddt(Pn) += ddtPn_perp_conduction; } else { - ddt(Pn) += - (2. / 3) - * Div_a_Grad_perp_flows(kappa_n_perp, Tn, ef_cond_perp_xlow, ef_cond_perp_ylow); + ddtPn_perp_conduction = (2. / 3) + * Div_a_Grad_perp_flows(kappa_n_perp, Tn, ef_cond_perp_xlow, + ef_cond_perp_ylow); + ddt(Pn) += ddtPn_perp_conduction; } // The factor here is likely 3/2 as this is pure energy flow, but needs checking. @@ -758,19 +760,21 @@ void NeutralMixed::finally(const Options& state) { // Neutral momentum TRACE("Neutral momentum"); - ddt(NVn) = -AA + ddtNVn_par_advection = -AA * FV::Div_par_fvv( // Momentum flow - Nnlim, Vn, sound_speed) + Nnlim, Vn, sound_speed); + ddt(NVn) += ddtNVn_par_advection; - - Grad_par(Pn); // Pressure gradient + ddtNVn_pressure_gradient = -Grad_par(Pn); // Pressure gradient + ddt(NVn) += ddtNVn_pressure_gradient; // Perpendicular advection of momentum if (nonorthogonal_operators) { - ddt(NVn) += - Div_a_Grad_perp_nonorthog(DnnNVn, logPnlim, mf_adv_perp_xlow, mf_adv_perp_ylow); + ddtNVn_perp_advection = Div_a_Grad_perp_nonorthog(DnnNVn, logPnlim, mf_adv_perp_xlow, mf_adv_perp_ylow); + ddt(NVn) += ddtNVn_perp_advection; } else { - ddt(NVn) += - Div_a_Grad_perp_flows(DnnNVn, logPnlim, mf_adv_perp_xlow, mf_adv_perp_ylow); + ddtNVn_perp_advection = Div_a_Grad_perp_flows(DnnNVn, logPnlim, mf_adv_perp_xlow, mf_adv_perp_ylow); + ddt(NVn) += ddtNVn_perp_advection; } if (neutral_viscosity) { @@ -782,131 +786,23 @@ void NeutralMixed::finally(const Options& state) { // Transport Processes in Gases", 1972 // eta_n = (2. / 5) * kappa_n; - Field3D viscosity_source = Div_par_K_Grad_par_mod( // Parallel viscosity + par_viscosity_source = Div_par_K_Grad_par_mod( // Parallel viscosity eta_n_par, Vn, mf_visc_par_ylow, false) // No viscosity through target boundary ; // Perpendicular viscosity if (nonorthogonal_operators) { - viscosity_source += Div_a_Grad_perp_nonorthog(eta_n_perp, Vn, mf_visc_perp_xlow, + perp_viscosity_source = Div_a_Grad_perp_nonorthog(eta_n_perp, Vn, mf_visc_perp_xlow, mf_visc_perp_ylow); } else { - viscosity_source += - Div_a_Grad_perp_flows(eta_n_perp, Vn, mf_visc_perp_xlow, mf_visc_perp_ylow); + perp_viscosity_source = Div_a_Grad_perp_flows(eta_n_perp, Vn, mf_visc_perp_xlow, mf_visc_perp_ylow); } - ddt(NVn) += viscosity_source; - ddt(Pn) += -(2. / 3) * Vn * viscosity_source; + ddt(NVn) += par_viscosity_source + perp_viscosity_source; + ddt(Pn) += -(2. / 3) * Vn * (par_viscosity_source + perp_viscosity_source); } Snv = momentum_source; - if (localstate.isSet("momentum_source")) { - Snv += get(localstate["momentum_source"]); - } - ddt(NVn) += Snv; - - } else { - ddt(NVn) = 0; - Snv = 0; - } - - // Add the contribution of ion perp velocity (i.e. anomalous transport) - // See eq 20 and 21 by Horsten et al., (2017) - if (perp_ion_coupling) { - - const Options& allspecies = state["species"]; - - for (auto& kv : allspecies.getChildren()) { - // NOTE:: This is only true for d+ ions. How do we generalize? - // How do we include the perpendicular ion velocity from other drifts? - - const Options& species = kv.second; - - if ((kv.first == "e") or !species.isSet("charge") - or (fabs(get(species["charge"])) < 1e-5)) { - continue; // Skip electrons and non-charged ions - } - - // sources/sinks due to anomalous transport - if (species.isSet("anomalous_D")) { - const Field2D anomalous_D = get(species["anomalous_D"]); - - const Field3D Ni = get(species["density"]); - Field2D Ni2D = DC(Ni); - - // Apply Neumann Y boundary condition, so no additional flux into boundary - // Note: Not setting radial (X) boundaries since those set radial fluxes - for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { - Ni2D(r.ind, mesh->ystart - 1) = Ni2D(r.ind, mesh->ystart); - } - for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { - Ni2D(r.ind, mesh->yend + 1) = Ni2D(r.ind, mesh->yend); - } - - ddt(Nn) += - Div_a_Grad_perp_upwind(Nn * anomalous_D / softFloor(Ni, density_floor), Ni2D); - // NOTE: Here, we used Nn as is done in UEDGE but it supposted to be the - // equilibrium value of Nn. - - ddt(Pn) += (5. / 3) - * Div_a_Grad_perp_upwind( - Pn * anomalous_D / softFloor(Ni, density_floor), Ni2D); - - if (evolve_momentum) { - ddt(NVn) += Div_a_Grad_perp_upwind( - NVn * anomalous_D / softFloor(Ni, density_floor), Ni2D); - } - } - } - } - - // If N < density_floor then NV and NV_solver may differ - // -> Add term to force NV_solver towards NV - // Note: This correction is calculated in transform() - ddt(NVn) += NVn_err; - - // Ste time derivatives to zero - if (zero_timederivs) { - - Field3D zero{0.0}; - zero.splitParallelSlices(); - zero.yup() = 0.0; - zero.ydown() = 0.0; - - ddt(Nn) = zero; - ddt(Pn) = zero; - if (evolve_momentum) { - ddt(NVn) = zero; - } - return; - } - - // Scale time derivatives - if (state.isSet("scale_timederivs")) { - Field3D scale_timederivs = get(state["scale_timederivs"]); - ddt(Nn) *= scale_timederivs; - ddt(Pn) *= scale_timederivs; - ddt(NVn) *= scale_timederivs; - } - - if (freeze_low_density) { - // Apply a factor to time derivatives in low density regions. - // Keep the sources and sinks, so that temperature and flow - // equilibriates with the plasma through collisions. - - Field3D Nn_s, Pn_s, NVn_s; - if (localstate.isSet("density_source")) { - Nn_s = get(localstate["density_source"]); - } else { - Nn_s = 0.0; - } - if (localstate.isSet("energy_source")) { - Pn_s = (2. / 3) * get(localstate["energy_source"]); - } else { - Pn_s = 0.0; - } - if (localstate.isSet("momentum_source")) { - NVn_s = get(localstate["momentum_source"]); } else { NVn_s = 0.0; } From 161f8bff3dd989b5f44371cb7bc3f870a65bc7be Mon Sep 17 00:00:00 2001 From: Lin Shih Date: Wed, 1 Apr 2026 11:41:16 +0100 Subject: [PATCH 03/12] add terms --- include/neutral_mixed.hxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index 2ca455452..aac292f43 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -108,14 +108,14 @@ private: bool output_ddt; ///< Save time derivatives? bool diagnose; ///< Save additional diagnostics? - //Physics terms + //Physic terms Field3D ddtN_par_advection; Field3D ddtN_perp_diffusion; Field3D ddtPn_par_advection; Field3D ddtPn_work_done; Field3D ddtPn_perp_advection; - Field3D ddtPn_par_conduction;:w + Field3D ddtPn_par_conduction; Field3D ddtPn_perp_conduction; Field3D ddtNVn_par_advection; From 60bdbc622b0510dd256367104664abf52a86c53a Mon Sep 17 00:00:00 2001 From: Lin Shih Date: Wed, 1 Apr 2026 16:29:05 +0100 Subject: [PATCH 04/12] add all diagnostic for each evolution term --- src/neutral_mixed.cxx | 183 +++++++++++++++++++++++++++++++++++++++++- 1 file changed, 182 insertions(+), 1 deletion(-) diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 99301945f..2a9569ee1 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -702,7 +702,7 @@ void NeutralMixed::finally(const Options& state) { ddtPn_par_advection = -(5. / 3) * FV::Div_par_mod( // Parallel advection - Pn, Vn, sound_speed, ef_adv_par_ylow) + Pn, Vn, sound_speed, ef_adv_par_ylow); ddtPn_work_done = (2. / 3) * Vn * Grad_par(Pn); // Work done ddt(Pn) = ddtPn_par_advection + ddtPn_work_done; @@ -803,6 +803,113 @@ void NeutralMixed::finally(const Options& state) { ddt(Pn) += -(2. / 3) * Vn * (par_viscosity_source + perp_viscosity_source); } Snv = momentum_source; + if (localstate.isSet("momentum_source")) { + Snv += get(localstate["momentum_source"]); + } + ddt(NVn) += Snv; + + } else { + ddt(NVn) = 0; + Snv = 0; + } + + // Add the contribution of ion perp velocity (i.e. anomalous transport) + // See eq 20 and 21 by Horsten et al., (2017) + if (perp_ion_coupling) { + + const Options& allspecies = state["species"]; + + for (auto& kv : allspecies.getChildren()) { + // NOTE:: This is only true for d+ ions. How do we generalize? + // How do we include the perpendicular ion velocity from other drifts? + + const Options& species = kv.second; + + if ((kv.first == "e") or !species.isSet("charge") + or (fabs(get(species["charge"])) < 1e-5)) { + continue; // Skip electrons and non-charged ions + } + + // sources/sinks due to anomalous transport + if (species.isSet("anomalous_D")) { + const Field2D anomalous_D = get(species["anomalous_D"]); + + const Field3D Ni = get(species["density"]); + Field2D Ni2D = DC(Ni); + + // Apply Neumann Y boundary condition, so no additional flux into boundary + // Note: Not setting radial (X) boundaries since those set radial fluxes + for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { + Ni2D(r.ind, mesh->ystart - 1) = Ni2D(r.ind, mesh->ystart); + } + for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { + Ni2D(r.ind, mesh->yend + 1) = Ni2D(r.ind, mesh->yend); + } + + ddt(Nn) += + Div_a_Grad_perp_upwind(Nn * anomalous_D / softFloor(Ni, density_floor), Ni2D); + // NOTE: Here, we used Nn as is done in UEDGE but it supposted to be the + // equilibrium value of Nn. + + ddt(Pn) += (5. / 3) + * Div_a_Grad_perp_upwind( + Pn * anomalous_D / softFloor(Ni, density_floor), Ni2D); + + if (evolve_momentum) { + ddt(NVn) += Div_a_Grad_perp_upwind( + NVn * anomalous_D / softFloor(Ni, density_floor), Ni2D); + } + } + } + } + + // If N < density_floor then NV and NV_solver may differ + // -> Add term to force NV_solver towards NV + // Note: This correction is calculated in transform() + ddt(NVn) += NVn_err; + + // Ste time derivatives to zero + if (zero_timederivs) { + + Field3D zero{0.0}; + zero.splitParallelSlices(); + zero.yup() = 0.0; + zero.ydown() = 0.0; + + ddt(Nn) = zero; + ddt(Pn) = zero; + if (evolve_momentum) { + ddt(NVn) = zero; + } + return; + } + + // Scale time derivatives + if (state.isSet("scale_timederivs")) { + Field3D scale_timederivs = get(state["scale_timederivs"]); + ddt(Nn) *= scale_timederivs; + ddt(Pn) *= scale_timederivs; + ddt(NVn) *= scale_timederivs; + } + + if (freeze_low_density) { + // Apply a factor to time derivatives in low density regions. + // Keep the sources and rescues, so that temperature and flow + // equilibriates with the plasma through collisions. + + Field3D Nn_s, Pn_s, NVn_s; + if (localstate.isSet("density_source")) { + Nn_s = get(localstate["density_source"]); + } else { + Nn_s = 0.0; + } + if (localstate.isSet("energy_source")) { + Pn_s = (2. / 3) * get(localstate["energy_source"]); + } else { + Pn_s = 0.0; + } + if (localstate.isSet("momentum_source")) { + NVn_s = get(localstate["momentum_source"]); } else { NVn_s = 0.0; } @@ -892,16 +999,90 @@ void NeutralMixed::outputVars(Options& state) { {"conversion", Nnorm * Omega_ci}, {"long_name", std::string("Rate of change of ") + name + " number density"}, {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("ddtN") + name + std::string("_par_advection")], ddtN_par_advection, + {{"time_dimension", "t"}, + {"units", "m^-3 s^-1"}, + {"conversion", Nnorm * Omega_ci}, + {"long_name", name + std::string(" density parallel advection")}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("ddtN") + name + std::string("_perp_diffusion")], ddtN_perp_diffusion, + {{"time_dimension", "t"}, + {"units", "m^-3 s^-1"}, + {"conversion", Nnorm * Omega_ci}, + {"long_name", name + std::string(" density perpendicular diffusion")}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("ddt(P") + name + std::string(")")], ddt(Pn), {{"time_dimension", "t"}, {"units", "Pa s^-1"}, {"conversion", Pnorm * Omega_ci}, {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("ddtP") + name + std::string("_par_advection")], ddtPn_par_advection, + {{"time_dimension", "t"}, + {"units", "Pa s^-1"}, + {"conversion", Pnorm * Omega_ci}, + {"long_name", name + std::string(" pressure parallel advection")}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("ddtP") + name + std::string("_work_done")], ddtPn_work_done, + {{"time_dimension", "t"}, + {"units", "Pa s^-1"}, + {"conversion", Pnorm * Omega_ci}, + {"long_name", name + std::string(" pressure work done")}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("ddtP") + name + std::string("_perp_advection")], ddtPn_perp_advection, + {{"time_dimension", "t"}, + {"units", "Pa s^-1"}, + {"conversion", Pnorm * Omega_ci}, + {"long_name", name + std::string(" pressure perpendicular advection")}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("ddtP") + name + std::string("_par_conduction")], ddtPn_par_conduction, + {{"time_dimension", "t"}, + {"units", "Pa s^-1"}, + {"conversion", Pnorm * Omega_ci}, + {"long_name", name + std::string(" pressure parallel conduction")}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("ddtP") + name + std::string("_perp_conduction")], ddtPn_perp_conduction, + {{"time_dimension", "t"}, + {"units", "Pa s^-1"}, + {"conversion", Pnorm * Omega_ci}, + {"long_name", name + std::string(" pressure perpendicular conduction")}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("ddt(NV") + name + std::string(")")], ddt(NVn), {{"time_dimension", "t"}, {"units", "kg m^-2 s^-2"}, {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("ddtNV") + name + std::string("_par_advection")], ddtNVn_par_advection, + {{"time_dimension", "t"}, + {"units", "kg m^-2 s^-2"}, + {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, + {"long_name", name + std::string(" momentum parallel advection")}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("ddtNV") + name + std::string("_pressure_gradient")], ddtNVn_pressure_gradient, + {{"time_dimension", "t"}, + {"units", "kg m^-2 s^-2"}, + {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, + {"long_name", name + std::string(" momentum pressure gradient")}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("ddtNV") + name + std::string("_perp_advection")], ddtNVn_perp_advection, + {{"time_dimension", "t"}, + {"units", "kg m^-2 s^-2"}, + {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, + {"long_name", name + std::string(" momentum perpendicular advection")}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("NV") + name + std::string("_par_visc_source")], par_viscosity_source, + {{"time_dimension", "t"}, + {"units", "kg m^-2 s^-2"}, + {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, + {"long_name", name + std::string(" momentum parallel viscosity source")}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("NV") + name + std::string("_perp_visc_source")], perp_viscosity_source, + {{"time_dimension", "t"}, + {"units", "kg m^-2 s^-2"}, + {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, + {"long_name", name + std::string(" momentum perpendicular viscosity source")}, + {"source", "neutral_mixed"}}); } if (diagnose) { set_with_attrs(state["debug"], debug, From 873c934dd9d5cb248270f93061847bd85383b76d Mon Sep 17 00:00:00 2001 From: Lin Shih Date: Thu, 9 Apr 2026 11:10:29 +0100 Subject: [PATCH 05/12] fix ddtNVn_par_advection += to = --- src/neutral_mixed.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 2a9569ee1..b81bd287f 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -763,7 +763,7 @@ void NeutralMixed::finally(const Options& state) { ddtNVn_par_advection = -AA * FV::Div_par_fvv( // Momentum flow Nnlim, Vn, sound_speed); - ddt(NVn) += ddtNVn_par_advection; + ddt(NVn) = ddtNVn_par_advection; ddtNVn_pressure_gradient = -Grad_par(Pn); // Pressure gradient ddt(NVn) += ddtNVn_pressure_gradient; From 0210aa3ef9bb3a9eb67c101c1242cf0a212a8862 Mon Sep 17 00:00:00 2001 From: Lin Shih Date: Fri, 10 Apr 2026 13:09:06 +0100 Subject: [PATCH 06/12] add viscosity term --- include/neutral_mixed.hxx | 2 ++ src/neutral_mixed.cxx | 6 ++++-- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index aac292f43..4ffe38f90 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -123,6 +123,8 @@ private: Field3D ddtNVn_perp_advection; Field3D par_viscosity_source; Field3D perp_viscosity_source; + Field3D ddtNVn_viscosity; // par_viscosity_source + perp_viscosity_source + Field3D ddtPn_viscosity; // Flow diagnostics Field3D pf_adv_perp_xlow, pf_adv_perp_ylow, pf_adv_par_ylow; diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index b81bd287f..24e24cd2d 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -799,8 +799,10 @@ void NeutralMixed::finally(const Options& state) { perp_viscosity_source = Div_a_Grad_perp_flows(eta_n_perp, Vn, mf_visc_perp_xlow, mf_visc_perp_ylow); } - ddt(NVn) += par_viscosity_source + perp_viscosity_source; - ddt(Pn) += -(2. / 3) * Vn * (par_viscosity_source + perp_viscosity_source); + ddtNVn_viscosity = par_viscosity_source + perp_viscosity_source; + ddtPn_viscosity = -(2. / 3) * Vn * (par_viscosity_source + perp_viscosity_source); + ddt(NVn) += ddtNVn_viscosity; + ddt(Pn) += ddtPn_viscosity; } Snv = momentum_source; if (localstate.isSet("momentum_source")) { From d5fc9f3f5f0a75b8de0f35aaafe541114a241951 Mon Sep 17 00:00:00 2001 From: sun51027 <60597328+sun51027@users.noreply.github.com> Date: Fri, 10 Apr 2026 12:21:33 +0000 Subject: [PATCH 07/12] [bot] Apply format changes --- include/neutral_mixed.hxx | 4 +- src/neutral_mixed.cxx | 176 ++++++++++++++++++++++---------------- 2 files changed, 102 insertions(+), 78 deletions(-) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index 4ffe38f90..7af7206a8 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -108,7 +108,7 @@ private: bool output_ddt; ///< Save time derivatives? bool diagnose; ///< Save additional diagnostics? - //Physic terms + // Physic terms Field3D ddtN_par_advection; Field3D ddtN_perp_diffusion; @@ -124,7 +124,7 @@ private: Field3D par_viscosity_source; Field3D perp_viscosity_source; Field3D ddtNVn_viscosity; // par_viscosity_source + perp_viscosity_source - Field3D ddtPn_viscosity; + Field3D ddtPn_viscosity; // Flow diagnostics Field3D pf_adv_perp_xlow, pf_adv_perp_ylow, pf_adv_par_ylow; diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 24e24cd2d..210ac670f 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -677,16 +677,17 @@ void NeutralMixed::finally(const Options& state) { ///////////////////////////////////////////////////// // Neutral density TRACE("Neutral density"); - ddtN_par_advection = -FV::Div_par_mod(Nn, Vn, sound_speed, - pf_adv_par_ylow); + ddtN_par_advection = -FV::Div_par_mod(Nn, Vn, sound_speed, pf_adv_par_ylow); ddt(Nn) = ddtN_par_advection; // Parallel advection // Perpendicular diffusion if (nonorthogonal_operators) { - ddtN_perp_diffusion = Div_a_Grad_perp_nonorthog(DnnNn, logPnlim, pf_adv_perp_xlow, pf_adv_perp_ylow); + ddtN_perp_diffusion = + Div_a_Grad_perp_nonorthog(DnnNn, logPnlim, pf_adv_perp_xlow, pf_adv_perp_ylow); ddt(Nn) += ddtN_perp_diffusion; } else { - ddtN_perp_diffusion = Div_a_Grad_perp_flows(DnnNn, logPnlim, pf_adv_perp_xlow, pf_adv_perp_ylow); + ddtN_perp_diffusion = + Div_a_Grad_perp_flows(DnnNn, logPnlim, pf_adv_perp_xlow, pf_adv_perp_ylow); ddt(Nn) += ddtN_perp_diffusion; } @@ -701,18 +702,20 @@ void NeutralMixed::finally(const Options& state) { TRACE("Neutral pressure"); ddtPn_par_advection = -(5. / 3) - * FV::Div_par_mod( // Parallel advection - Pn, Vn, sound_speed, ef_adv_par_ylow); + * FV::Div_par_mod( // Parallel advection + Pn, Vn, sound_speed, ef_adv_par_ylow); ddtPn_work_done = (2. / 3) * Vn * Grad_par(Pn); // Work done ddt(Pn) = ddtPn_par_advection + ddtPn_work_done; // Perpendicular advection of pressure if (nonorthogonal_operators) { - ddtPn_perp_advection = (5. / 3) + ddtPn_perp_advection = + (5. / 3) * Div_a_Grad_perp_nonorthog(DnnPn, logPnlim, ef_adv_perp_xlow, ef_adv_perp_ylow); ddt(Pn) += ddtPn_perp_advection; } else { - ddtPn_perp_advection = (5. / 3) + ddtPn_perp_advection = + (5. / 3) * Div_a_Grad_perp_flows(DnnPn, logPnlim, ef_adv_perp_xlow, ef_adv_perp_ylow); ddt(Pn) += ddtPn_perp_advection; } @@ -723,22 +726,23 @@ void NeutralMixed::finally(const Options& state) { ef_adv_perp_ylow *= 5. / 2; if (neutral_conduction) { - ddtPn_par_conduction = (2. / 3) - * Div_par_K_Grad_par_mod(kappa_n_par, Tn, // Parallel conduction - ef_cond_par_ylow, - false); // No conduction through target boundary + ddtPn_par_conduction = + (2. / 3) + * Div_par_K_Grad_par_mod(kappa_n_par, Tn, // Parallel conduction + ef_cond_par_ylow, + false); // No conduction through target boundary ddt(Pn) += ddtPn_par_conduction; // Perpendicular conduction if (nonorthogonal_operators) { ddtPn_perp_conduction = (2. / 3) - * Div_a_Grad_perp_nonorthog(kappa_n_perp, Tn, ef_cond_perp_xlow, - ef_cond_perp_ylow); + * Div_a_Grad_perp_nonorthog( + kappa_n_perp, Tn, ef_cond_perp_xlow, ef_cond_perp_ylow); ddt(Pn) += ddtPn_perp_conduction; } else { - ddtPn_perp_conduction = (2. / 3) - * Div_a_Grad_perp_flows(kappa_n_perp, Tn, ef_cond_perp_xlow, - ef_cond_perp_ylow); + ddtPn_perp_conduction = + (2. / 3) + * Div_a_Grad_perp_flows(kappa_n_perp, Tn, ef_cond_perp_xlow, ef_cond_perp_ylow); ddt(Pn) += ddtPn_perp_conduction; } @@ -761,8 +765,8 @@ void NeutralMixed::finally(const Options& state) { TRACE("Neutral momentum"); ddtNVn_par_advection = -AA - * FV::Div_par_fvv( // Momentum flow - Nnlim, Vn, sound_speed); + * FV::Div_par_fvv( // Momentum flow + Nnlim, Vn, sound_speed); ddt(NVn) = ddtNVn_par_advection; ddtNVn_pressure_gradient = -Grad_par(Pn); // Pressure gradient @@ -770,10 +774,12 @@ void NeutralMixed::finally(const Options& state) { // Perpendicular advection of momentum if (nonorthogonal_operators) { - ddtNVn_perp_advection = Div_a_Grad_perp_nonorthog(DnnNVn, logPnlim, mf_adv_perp_xlow, mf_adv_perp_ylow); + ddtNVn_perp_advection = + Div_a_Grad_perp_nonorthog(DnnNVn, logPnlim, mf_adv_perp_xlow, mf_adv_perp_ylow); ddt(NVn) += ddtNVn_perp_advection; } else { - ddtNVn_perp_advection = Div_a_Grad_perp_flows(DnnNVn, logPnlim, mf_adv_perp_xlow, mf_adv_perp_ylow); + ddtNVn_perp_advection = + Div_a_Grad_perp_flows(DnnNVn, logPnlim, mf_adv_perp_xlow, mf_adv_perp_ylow); ddt(NVn) += ddtNVn_perp_advection; } @@ -793,16 +799,17 @@ void NeutralMixed::finally(const Options& state) { // Perpendicular viscosity if (nonorthogonal_operators) { - perp_viscosity_source = Div_a_Grad_perp_nonorthog(eta_n_perp, Vn, mf_visc_perp_xlow, - mf_visc_perp_ylow); + perp_viscosity_source = Div_a_Grad_perp_nonorthog( + eta_n_perp, Vn, mf_visc_perp_xlow, mf_visc_perp_ylow); } else { - perp_viscosity_source = Div_a_Grad_perp_flows(eta_n_perp, Vn, mf_visc_perp_xlow, mf_visc_perp_ylow); + perp_viscosity_source = + Div_a_Grad_perp_flows(eta_n_perp, Vn, mf_visc_perp_xlow, mf_visc_perp_ylow); } - ddtNVn_viscosity = par_viscosity_source + perp_viscosity_source; - ddtPn_viscosity = -(2. / 3) * Vn * (par_viscosity_source + perp_viscosity_source); + ddtNVn_viscosity = par_viscosity_source + perp_viscosity_source; + ddtPn_viscosity = -(2. / 3) * Vn * (par_viscosity_source + perp_viscosity_source); ddt(NVn) += ddtNVn_viscosity; - ddt(Pn) += ddtPn_viscosity; + ddt(Pn) += ddtPn_viscosity; } Snv = momentum_source; if (localstate.isSet("momentum_source")) { @@ -1001,13 +1008,15 @@ void NeutralMixed::outputVars(Options& state) { {"conversion", Nnorm * Omega_ci}, {"long_name", std::string("Rate of change of ") + name + " number density"}, {"source", "neutral_mixed"}}); - set_with_attrs(state[std::string("ddtN") + name + std::string("_par_advection")], ddtN_par_advection, + set_with_attrs(state[std::string("ddtN") + name + std::string("_par_advection")], + ddtN_par_advection, {{"time_dimension", "t"}, {"units", "m^-3 s^-1"}, {"conversion", Nnorm * Omega_ci}, {"long_name", name + std::string(" density parallel advection")}, {"source", "neutral_mixed"}}); - set_with_attrs(state[std::string("ddtN") + name + std::string("_perp_diffusion")], ddtN_perp_diffusion, + set_with_attrs(state[std::string("ddtN") + name + std::string("_perp_diffusion")], + ddtN_perp_diffusion, {{"time_dimension", "t"}, {"units", "m^-3 s^-1"}, {"conversion", Nnorm * Omega_ci}, @@ -1019,72 +1028,87 @@ void NeutralMixed::outputVars(Options& state) { {"units", "Pa s^-1"}, {"conversion", Pnorm * Omega_ci}, {"source", "neutral_mixed"}}); - set_with_attrs(state[std::string("ddtP") + name + std::string("_par_advection")], ddtPn_par_advection, + set_with_attrs(state[std::string("ddtP") + name + std::string("_par_advection")], + ddtPn_par_advection, {{"time_dimension", "t"}, {"units", "Pa s^-1"}, {"conversion", Pnorm * Omega_ci}, {"long_name", name + std::string(" pressure parallel advection")}, {"source", "neutral_mixed"}}); - set_with_attrs(state[std::string("ddtP") + name + std::string("_work_done")], ddtPn_work_done, + set_with_attrs(state[std::string("ddtP") + name + std::string("_work_done")], + ddtPn_work_done, {{"time_dimension", "t"}, {"units", "Pa s^-1"}, {"conversion", Pnorm * Omega_ci}, {"long_name", name + std::string(" pressure work done")}, {"source", "neutral_mixed"}}); - set_with_attrs(state[std::string("ddtP") + name + std::string("_perp_advection")], ddtPn_perp_advection, - {{"time_dimension", "t"}, - {"units", "Pa s^-1"}, - {"conversion", Pnorm * Omega_ci}, - {"long_name", name + std::string(" pressure perpendicular advection")}, - {"source", "neutral_mixed"}}); - set_with_attrs(state[std::string("ddtP") + name + std::string("_par_conduction")], ddtPn_par_conduction, + set_with_attrs( + state[std::string("ddtP") + name + std::string("_perp_advection")], + ddtPn_perp_advection, + {{"time_dimension", "t"}, + {"units", "Pa s^-1"}, + {"conversion", Pnorm * Omega_ci}, + {"long_name", name + std::string(" pressure perpendicular advection")}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("ddtP") + name + std::string("_par_conduction")], + ddtPn_par_conduction, {{"time_dimension", "t"}, {"units", "Pa s^-1"}, {"conversion", Pnorm * Omega_ci}, {"long_name", name + std::string(" pressure parallel conduction")}, {"source", "neutral_mixed"}}); - set_with_attrs(state[std::string("ddtP") + name + std::string("_perp_conduction")], ddtPn_perp_conduction, - {{"time_dimension", "t"}, - {"units", "Pa s^-1"}, - {"conversion", Pnorm * Omega_ci}, - {"long_name", name + std::string(" pressure perpendicular conduction")}, - {"source", "neutral_mixed"}}); + set_with_attrs( + state[std::string("ddtP") + name + std::string("_perp_conduction")], + ddtPn_perp_conduction, + {{"time_dimension", "t"}, + {"units", "Pa s^-1"}, + {"conversion", Pnorm * Omega_ci}, + {"long_name", name + std::string(" pressure perpendicular conduction")}, + {"source", "neutral_mixed"}}); set_with_attrs(state[std::string("ddt(NV") + name + std::string(")")], ddt(NVn), {{"time_dimension", "t"}, {"units", "kg m^-2 s^-2"}, {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, {"source", "neutral_mixed"}}); - set_with_attrs(state[std::string("ddtNV") + name + std::string("_par_advection")], ddtNVn_par_advection, + set_with_attrs(state[std::string("ddtNV") + name + std::string("_par_advection")], + ddtNVn_par_advection, {{"time_dimension", "t"}, {"units", "kg m^-2 s^-2"}, {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, {"long_name", name + std::string(" momentum parallel advection")}, {"source", "neutral_mixed"}}); - set_with_attrs(state[std::string("ddtNV") + name + std::string("_pressure_gradient")], ddtNVn_pressure_gradient, + set_with_attrs(state[std::string("ddtNV") + name + std::string("_pressure_gradient")], + ddtNVn_pressure_gradient, {{"time_dimension", "t"}, {"units", "kg m^-2 s^-2"}, {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, {"long_name", name + std::string(" momentum pressure gradient")}, {"source", "neutral_mixed"}}); - set_with_attrs(state[std::string("ddtNV") + name + std::string("_perp_advection")], ddtNVn_perp_advection, - {{"time_dimension", "t"}, - {"units", "kg m^-2 s^-2"}, - {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, - {"long_name", name + std::string(" momentum perpendicular advection")}, - {"source", "neutral_mixed"}}); - set_with_attrs(state[std::string("NV") + name + std::string("_par_visc_source")], par_viscosity_source, - {{"time_dimension", "t"}, - {"units", "kg m^-2 s^-2"}, - {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, - {"long_name", name + std::string(" momentum parallel viscosity source")}, - {"source", "neutral_mixed"}}); - set_with_attrs(state[std::string("NV") + name + std::string("_perp_visc_source")], perp_viscosity_source, - {{"time_dimension", "t"}, - {"units", "kg m^-2 s^-2"}, - {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, - {"long_name", name + std::string(" momentum perpendicular viscosity source")}, - {"source", "neutral_mixed"}}); + set_with_attrs( + state[std::string("ddtNV") + name + std::string("_perp_advection")], + ddtNVn_perp_advection, + {{"time_dimension", "t"}, + {"units", "kg m^-2 s^-2"}, + {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, + {"long_name", name + std::string(" momentum perpendicular advection")}, + {"source", "neutral_mixed"}}); + set_with_attrs( + state[std::string("NV") + name + std::string("_par_visc_source")], + par_viscosity_source, + {{"time_dimension", "t"}, + {"units", "kg m^-2 s^-2"}, + {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, + {"long_name", name + std::string(" momentum parallel viscosity source")}, + {"source", "neutral_mixed"}}); + set_with_attrs( + state[std::string("NV") + name + std::string("_perp_visc_source")], + perp_viscosity_source, + {{"time_dimension", "t"}, + {"units", "kg m^-2 s^-2"}, + {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, + {"long_name", name + std::string(" momentum perpendicular viscosity source")}, + {"source", "neutral_mixed"}}); } if (diagnose) { set_with_attrs(state["debug"], debug, @@ -1247,19 +1271,19 @@ void NeutralMixed::outputVars(Options& state) { {"source", "neutral_mixed"}}); // Lin add set_with_attrs(state[std::string("V") + name + std::string("th_pf")], Vnth_pf, - {{"time_dimension", "t"}, - {"units", "m / s"}, - {"conversion", Cs0}, - {"standard_name", "velocity"}, - {"long_name", name + " thermal velocity of neutrals for advection"}, - {"source", "neutral_mixed"}}); + {{"time_dimension", "t"}, + {"units", "m / s"}, + {"conversion", Cs0}, + {"standard_name", "velocity"}, + {"long_name", name + " thermal velocity of neutrals for advection"}, + {"source", "neutral_mixed"}}); set_with_attrs(state[std::string("V") + name + std::string("th_hf")], Vnth_hf, - {{"time_dimension", "t"}, - {"units", "m / s"}, - {"conversion", Cs0}, - {"standard_name", "velocity"}, - {"long_name", name + " thermal velocity of neutrals for heat flux"}, - {"source", "neutral_mixed"}}); + {{"time_dimension", "t"}, + {"units", "m / s"}, + {"conversion", Cs0}, + {"standard_name", "velocity"}, + {"long_name", name + " thermal velocity of neutrals for heat flux"}, + {"source", "neutral_mixed"}}); /////////////////////////////////////////////////// // Parallel flow diagnostics From 96fbef608692007d9f50d23e06b836605f332b14 Mon Sep 17 00:00:00 2001 From: sun51027 <60597328+sun51027@users.noreply.github.com> Date: Fri, 10 Apr 2026 12:21:33 +0000 Subject: [PATCH 08/12] [bot] Add last format changes commit to ignore file --- .git-blame-ignore-revs | 1 + 1 file changed, 1 insertion(+) diff --git a/.git-blame-ignore-revs b/.git-blame-ignore-revs index 5ea4e3c5d..92411b684 100644 --- a/.git-blame-ignore-revs +++ b/.git-blame-ignore-revs @@ -2,3 +2,4 @@ f5de620399e88f139e236e40675facfd0cfd6b13 07616298519738fad0c1bb259e4d021cd6f51d58 37efa12abe0961add8715b797d3926664977d453 8e49d5933131ddcb2f97bb2ac7095a1be80ad7b7 +d5fc9f3f5f0a75b8de0f35aaafe541114a241951 From 75314e5e5ea04a332293a5b0075af322e8b9e9d9 Mon Sep 17 00:00:00 2001 From: Lin Shih Date: Fri, 10 Apr 2026 14:49:27 +0100 Subject: [PATCH 09/12] add viscosity and anomalous_D term --- include/neutral_mixed.hxx | 4 ++ src/neutral_mixed.cxx | 85 +++++++++++++++++++++++++-------------- 2 files changed, 58 insertions(+), 31 deletions(-) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index 7af7206a8..61779288c 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -126,6 +126,10 @@ private: Field3D ddtNVn_viscosity; // par_viscosity_source + perp_viscosity_source Field3D ddtPn_viscosity; + Field3D ddtN_anomalous_transport; + Field3D ddtNVn_anomalous_transport; + Field3D ddtPn_anomalous_transport; + // Flow diagnostics Field3D pf_adv_perp_xlow, pf_adv_perp_ylow, pf_adv_par_ylow; Field3D mf_adv_perp_xlow, mf_adv_perp_ylow, mf_adv_par_ylow; diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 210ac670f..b63392640 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -855,18 +855,17 @@ void NeutralMixed::finally(const Options& state) { Ni2D(r.ind, mesh->yend + 1) = Ni2D(r.ind, mesh->yend); } - ddt(Nn) += - Div_a_Grad_perp_upwind(Nn * anomalous_D / softFloor(Ni, density_floor), Ni2D); + ddtN_anomalous_transport = Div_a_Grad_perp_upwind(Nn * anomalous_D / softFloor(Ni, density_floor), Ni2D); + ddt(Nn) += ddtN_anomalous_transport; // NOTE: Here, we used Nn as is done in UEDGE but it supposted to be the // equilibrium value of Nn. - ddt(Pn) += (5. / 3) - * Div_a_Grad_perp_upwind( - Pn * anomalous_D / softFloor(Ni, density_floor), Ni2D); + ddtPn_anomalous_transport = (5. / 3) * Div_a_Grad_perp_upwind(Pn * anomalous_D / softFloor(Ni, density_floor), Ni2D); + ddt(Pn) += ddtPn_anomalous_transport; if (evolve_momentum) { - ddt(NVn) += Div_a_Grad_perp_upwind( - NVn * anomalous_D / softFloor(Ni, density_floor), Ni2D); + ddtNVn_anomalous_transport = Div_a_Grad_perp_upwind( NVn * anomalous_D / softFloor(Ni, density_floor), Ni2D); + ddt(NVn) += ddtNVn_anomalous_transport; } } } @@ -1085,30 +1084,54 @@ void NeutralMixed::outputVars(Options& state) { {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, {"long_name", name + std::string(" momentum pressure gradient")}, {"source", "neutral_mixed"}}); - set_with_attrs( - state[std::string("ddtNV") + name + std::string("_perp_advection")], - ddtNVn_perp_advection, - {{"time_dimension", "t"}, - {"units", "kg m^-2 s^-2"}, - {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, - {"long_name", name + std::string(" momentum perpendicular advection")}, - {"source", "neutral_mixed"}}); - set_with_attrs( - state[std::string("NV") + name + std::string("_par_visc_source")], - par_viscosity_source, - {{"time_dimension", "t"}, - {"units", "kg m^-2 s^-2"}, - {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, - {"long_name", name + std::string(" momentum parallel viscosity source")}, - {"source", "neutral_mixed"}}); - set_with_attrs( - state[std::string("NV") + name + std::string("_perp_visc_source")], - perp_viscosity_source, - {{"time_dimension", "t"}, - {"units", "kg m^-2 s^-2"}, - {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, - {"long_name", name + std::string(" momentum perpendicular viscosity source")}, - {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("ddtNV") + name + std::string("_perp_advection")], ddtNVn_perp_advection, + {{"time_dimension", "t"}, + {"units", "kg m^-2 s^-2"}, + {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, + {"long_name", name + std::string(" momentum perpendicular advection")}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("NV") + name + std::string("_par_visc_source")], par_viscosity_source, + {{"time_dimension", "t"}, + {"units", "kg m^-2 s^-2"}, + {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, + {"long_name", name + std::string(" momentum parallel viscosity source")}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("NV") + name + std::string("_perp_visc_source")], perp_viscosity_source, + {{"time_dimension", "t"}, + {"units", "kg m^-2 s^-2"}, + {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, + {"long_name", name + std::string(" momentum perpendicular viscosity source")}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("ddtNV") + name + std::string("_viscosity")], ddtNVn_viscosity, + {{"time_dimension", "t"}, + {"units", "kg m^-2 s^-2"}, + {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, + {"long_name", name + std::string(" momentum viscosity")}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("ddtPn") + name + std::string("_viscosity")], ddtPn_viscosity, + {{"time_dimension", "t"}, + {"units", "Pa s^-1"}, + {"conversion",Pnorm * Omega_ci}, + {"long_name", name + std::string(" pressure viscosity")}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("ddtN") + name + std::string("_anomalous_transport")], ddtN_anomalous_transport, + {{"time_dimension", "t"}, + {"units", "m^-3 s^-1"}, + {"conversion", Nnorm * Omega_ci}, + {"long_name", name + std::string(" density anomalous transport")}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("ddtNV") + name + std::string("_anomalous_transport")], ddtNVn_anomalous_transport, + {{"time_dimension", "t"}, + {"units", "kg m^-2 s^-2"}, + {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, + {"long_name", name + std::string(" momentum anomalous transport")}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("ddtPn") + name + std::string("_anomalous_transport")], ddtPn_anomalous_transport, + {{"time_dimension", "t"}, + {"units", "Pa s^-1"}, + {"conversion",Pnorm * Omega_ci}, + {"long_name", name + std::string(" pressure anomalous transport")}, + {"source", "neutral_mixed"}}); } if (diagnose) { set_with_attrs(state["debug"], debug, From 5eeba360b09f53a7b58fd3354bf349f36ab3d17f Mon Sep 17 00:00:00 2001 From: Lin Shih Date: Fri, 10 Apr 2026 20:27:42 +0100 Subject: [PATCH 10/12] initialise terms = 0.0 if boolean variable is required --- include/neutral_mixed.hxx | 43 ++++++++++++++++++++------------------- 1 file changed, 22 insertions(+), 21 deletions(-) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index 61779288c..1301e7452 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -108,27 +108,28 @@ private: bool output_ddt; ///< Save time derivatives? bool diagnose; ///< Save additional diagnostics? - // Physic terms - Field3D ddtN_par_advection; - Field3D ddtN_perp_diffusion; - - Field3D ddtPn_par_advection; - Field3D ddtPn_work_done; - Field3D ddtPn_perp_advection; - Field3D ddtPn_par_conduction; - Field3D ddtPn_perp_conduction; - - Field3D ddtNVn_par_advection; - Field3D ddtNVn_pressure_gradient; - Field3D ddtNVn_perp_advection; - Field3D par_viscosity_source; - Field3D perp_viscosity_source; - Field3D ddtNVn_viscosity; // par_viscosity_source + perp_viscosity_source - Field3D ddtPn_viscosity; - - Field3D ddtN_anomalous_transport; - Field3D ddtNVn_anomalous_transport; - Field3D ddtPn_anomalous_transport; + //Physic terms + Field3D ddtN_par_advection = 0.0; + Field3D ddtN_perp_diffusion = 0.0; + + Field3D ddtPn_par_advection = 0.0; + Field3D ddtPn_work_done = 0.0; + Field3D ddtPn_perp_advection = 0.0; + Field3D ddtPn_par_conduction = 0.0; + Field3D ddtPn_perp_conduction = 0.0; + + Field3D ddtNVn_par_advection = 0.0; + Field3D ddtNVn_pressure_gradient = 0.0; + Field3D ddtNVn_perp_advection = 0.0; + + Field3D par_viscosity_source = 0.0; + Field3D perp_viscosity_source = 0.0; + Field3D ddtNVn_viscosity = 0.0; // par_viscosity_source + perp_viscosity_source + Field3D ddtPn_viscosity = 0.0; + + Field3D ddtN_anomalous_transport = 0.0; + Field3D ddtNVn_anomalous_transport = 0.0; + Field3D ddtPn_anomalous_transport = 0.0; // Flow diagnostics Field3D pf_adv_perp_xlow, pf_adv_perp_ylow, pf_adv_par_ylow; From acf4625dd0940720ccfdc38af61493e19cde21ed Mon Sep 17 00:00:00 2001 From: sun51027 <60597328+sun51027@users.noreply.github.com> Date: Mon, 20 Apr 2026 14:23:22 +0000 Subject: [PATCH 11/12] [bot] Apply format changes --- include/neutral_mixed.hxx | 4 +- src/neutral_mixed.cxx | 109 ++++++++++++++++++++++---------------- 2 files changed, 66 insertions(+), 47 deletions(-) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index 1301e7452..ba8515d12 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -108,7 +108,7 @@ private: bool output_ddt; ///< Save time derivatives? bool diagnose; ///< Save additional diagnostics? - //Physic terms + // Physic terms Field3D ddtN_par_advection = 0.0; Field3D ddtN_perp_diffusion = 0.0; @@ -125,7 +125,7 @@ private: Field3D par_viscosity_source = 0.0; Field3D perp_viscosity_source = 0.0; Field3D ddtNVn_viscosity = 0.0; // par_viscosity_source + perp_viscosity_source - Field3D ddtPn_viscosity = 0.0; + Field3D ddtPn_viscosity = 0.0; Field3D ddtN_anomalous_transport = 0.0; Field3D ddtNVn_anomalous_transport = 0.0; diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index b63392640..6a8c97392 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -855,17 +855,22 @@ void NeutralMixed::finally(const Options& state) { Ni2D(r.ind, mesh->yend + 1) = Ni2D(r.ind, mesh->yend); } - ddtN_anomalous_transport = Div_a_Grad_perp_upwind(Nn * anomalous_D / softFloor(Ni, density_floor), Ni2D); - ddt(Nn) += ddtN_anomalous_transport; + ddtN_anomalous_transport = + Div_a_Grad_perp_upwind(Nn * anomalous_D / softFloor(Ni, density_floor), Ni2D); + ddt(Nn) += ddtN_anomalous_transport; // NOTE: Here, we used Nn as is done in UEDGE but it supposted to be the // equilibrium value of Nn. - ddtPn_anomalous_transport = (5. / 3) * Div_a_Grad_perp_upwind(Pn * anomalous_D / softFloor(Ni, density_floor), Ni2D); - ddt(Pn) += ddtPn_anomalous_transport; + ddtPn_anomalous_transport = + (5. / 3) + * Div_a_Grad_perp_upwind(Pn * anomalous_D / softFloor(Ni, density_floor), + Ni2D); + ddt(Pn) += ddtPn_anomalous_transport; if (evolve_momentum) { - ddtNVn_anomalous_transport = Div_a_Grad_perp_upwind( NVn * anomalous_D / softFloor(Ni, density_floor), Ni2D); - ddt(NVn) += ddtNVn_anomalous_transport; + ddtNVn_anomalous_transport = Div_a_Grad_perp_upwind( + NVn * anomalous_D / softFloor(Ni, density_floor), Ni2D); + ddt(NVn) += ddtNVn_anomalous_transport; } } } @@ -1084,54 +1089,68 @@ void NeutralMixed::outputVars(Options& state) { {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, {"long_name", name + std::string(" momentum pressure gradient")}, {"source", "neutral_mixed"}}); - set_with_attrs(state[std::string("ddtNV") + name + std::string("_perp_advection")], ddtNVn_perp_advection, - {{"time_dimension", "t"}, - {"units", "kg m^-2 s^-2"}, - {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, - {"long_name", name + std::string(" momentum perpendicular advection")}, - {"source", "neutral_mixed"}}); - set_with_attrs(state[std::string("NV") + name + std::string("_par_visc_source")], par_viscosity_source, - {{"time_dimension", "t"}, - {"units", "kg m^-2 s^-2"}, - {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, - {"long_name", name + std::string(" momentum parallel viscosity source")}, - {"source", "neutral_mixed"}}); - set_with_attrs(state[std::string("NV") + name + std::string("_perp_visc_source")], perp_viscosity_source, - {{"time_dimension", "t"}, - {"units", "kg m^-2 s^-2"}, - {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, - {"long_name", name + std::string(" momentum perpendicular viscosity source")}, - {"source", "neutral_mixed"}}); - set_with_attrs(state[std::string("ddtNV") + name + std::string("_viscosity")], ddtNVn_viscosity, + set_with_attrs( + state[std::string("ddtNV") + name + std::string("_perp_advection")], + ddtNVn_perp_advection, + {{"time_dimension", "t"}, + {"units", "kg m^-2 s^-2"}, + {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, + {"long_name", name + std::string(" momentum perpendicular advection")}, + {"source", "neutral_mixed"}}); + set_with_attrs( + state[std::string("NV") + name + std::string("_par_visc_source")], + par_viscosity_source, + {{"time_dimension", "t"}, + {"units", "kg m^-2 s^-2"}, + {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, + {"long_name", name + std::string(" momentum parallel viscosity source")}, + {"source", "neutral_mixed"}}); + set_with_attrs( + state[std::string("NV") + name + std::string("_perp_visc_source")], + perp_viscosity_source, + {{"time_dimension", "t"}, + {"units", "kg m^-2 s^-2"}, + {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, + {"long_name", name + std::string(" momentum perpendicular viscosity source")}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("ddtNV") + name + std::string("_viscosity")], + ddtNVn_viscosity, {{"time_dimension", "t"}, {"units", "kg m^-2 s^-2"}, {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, {"long_name", name + std::string(" momentum viscosity")}, {"source", "neutral_mixed"}}); - set_with_attrs(state[std::string("ddtPn") + name + std::string("_viscosity")], ddtPn_viscosity, + set_with_attrs(state[std::string("ddtPn") + name + std::string("_viscosity")], + ddtPn_viscosity, {{"time_dimension", "t"}, {"units", "Pa s^-1"}, - {"conversion",Pnorm * Omega_ci}, + {"conversion", Pnorm * Omega_ci}, {"long_name", name + std::string(" pressure viscosity")}, {"source", "neutral_mixed"}}); - set_with_attrs(state[std::string("ddtN") + name + std::string("_anomalous_transport")], ddtN_anomalous_transport, - {{"time_dimension", "t"}, - {"units", "m^-3 s^-1"}, - {"conversion", Nnorm * Omega_ci}, - {"long_name", name + std::string(" density anomalous transport")}, - {"source", "neutral_mixed"}}); - set_with_attrs(state[std::string("ddtNV") + name + std::string("_anomalous_transport")], ddtNVn_anomalous_transport, - {{"time_dimension", "t"}, - {"units", "kg m^-2 s^-2"}, - {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, - {"long_name", name + std::string(" momentum anomalous transport")}, - {"source", "neutral_mixed"}}); - set_with_attrs(state[std::string("ddtPn") + name + std::string("_anomalous_transport")], ddtPn_anomalous_transport, - {{"time_dimension", "t"}, - {"units", "Pa s^-1"}, - {"conversion",Pnorm * Omega_ci}, - {"long_name", name + std::string(" pressure anomalous transport")}, - {"source", "neutral_mixed"}}); + set_with_attrs( + state[std::string("ddtN") + name + std::string("_anomalous_transport")], + ddtN_anomalous_transport, + {{"time_dimension", "t"}, + {"units", "m^-3 s^-1"}, + {"conversion", Nnorm * Omega_ci}, + {"long_name", name + std::string(" density anomalous transport")}, + {"source", "neutral_mixed"}}); + set_with_attrs( + state[std::string("ddtNV") + name + std::string("_anomalous_transport")], + ddtNVn_anomalous_transport, + {{"time_dimension", "t"}, + {"units", "kg m^-2 s^-2"}, + {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, + {"long_name", name + std::string(" momentum anomalous transport")}, + {"source", "neutral_mixed"}}); + set_with_attrs( + state[std::string("ddtPn") + name + std::string("_anomalous_transport")], + ddtPn_anomalous_transport, + {{"time_dimension", "t"}, + {"units", "Pa s^-1"}, + {"conversion", Pnorm * Omega_ci}, + {"long_name", name + std::string(" pressure anomalous transport")}, + {"source", "neutral_mixed"}}); } if (diagnose) { set_with_attrs(state["debug"], debug, From 27c0bf8cc33d93d3caead19d3fb0af9d02181526 Mon Sep 17 00:00:00 2001 From: sun51027 <60597328+sun51027@users.noreply.github.com> Date: Mon, 20 Apr 2026 14:23:23 +0000 Subject: [PATCH 12/12] [bot] Add last format changes commit to ignore file --- .git-blame-ignore-revs | 1 + 1 file changed, 1 insertion(+) diff --git a/.git-blame-ignore-revs b/.git-blame-ignore-revs index 92411b684..93c365c6f 100644 --- a/.git-blame-ignore-revs +++ b/.git-blame-ignore-revs @@ -3,3 +3,4 @@ f5de620399e88f139e236e40675facfd0cfd6b13 37efa12abe0961add8715b797d3926664977d453 8e49d5933131ddcb2f97bb2ac7095a1be80ad7b7 d5fc9f3f5f0a75b8de0f35aaafe541114a241951 +acf4625dd0940720ccfdc38af61493e19cde21ed