Conversation
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #592 +/- ##
==========================================
- Coverage 57.73% 57.34% -0.39%
==========================================
Files 97 97
Lines 10013 10112 +99
Branches 1442 1475 +33
==========================================
+ Hits 5781 5799 +18
- Misses 3652 3725 +73
- Partials 580 588 +8 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
Re-applies the previously-proposed neutral_passive_momentum feature onto current master (the original branch had drifted ~600 commits behind). New options: - passive_momentum: evolve only the neutral density; momentum and pressure are passive (pressure is not added to the solver). - temperature_from: take Tn from another species (e.g. the ion, Tn=Ti). - evolve_pressure: defaults to !passive_momentum. - anomalous_transport_momentum (AnomalousDiffusion): toggle anomalous transport of parallel momentum. In passive_momentum mode NVn is set from an ion-equilibrium relation NVn = Nn_eq * Vi, with Nn_eq built from cx/iz/rec collision frequencies, and nu_cx/nu_iz/nu_rec/Nn_eq exposed as diagnostics. The preconditioner gains a density-only branch (I - gamma*Div_perp(Dnn*Grad_perp))*ddt(Nn) when pressure is not evolved. Pressure/momentum terms, sources, diagnostics, time-derivative scaling, freeze_low_density and CHECK blocks are gated on the new flags. Reading another species' temperature declares the required GuardedOptions read permission. Tested on a 3D MAST-U Super-X case: builds clean, runs to completion with finite fields, reproduces the original branch's passive behaviour (Tn tracks Ti to ~0.1%, 7 evolved fields, consistent neutral diagnostics).
29f1eb4 to
5aa5ae8
Compare
| solver->add(Nn, std::string("N") + name); | ||
| solver->add(Pn, std::string("P") + name); | ||
|
|
||
| evolve_momentum = options["evolve_momentum"] |
There was a problem hiding this comment.
It looks like when passive_momentum is enabled, you need to disable evolve_momentum manually or else it'll be evolved.
One solution could be to make evolve_momentum default to !passive_momentum, and raise a BoutException when both are enabled.
| if (!evolve_pressure && !temperature_from.empty()) { | ||
| Tn = get<Field3D>(localstate["temperature"]); | ||
| Pn = Nn * Tn; | ||
| } |
There was a problem hiding this comment.
What happens if evolve_pressure = False and temperature_from = ""?
| ///////////////////////////////////////////////////// | ||
| // Neutral density | ||
| TRACE("Neutral density"); | ||
| ddt(Nn) = -FV::Div_par_mod<ParLimiter>(Nn, Vn, sound_speed, |
There was a problem hiding this comment.
The density equation needs Vn. But at this point in the code, Vn is taken from the state, instead of the later recalculation from the ion velocity! This makes it inconsistent.
Usually both Vn and Nn are evolved by the solver in transform(). This means both the density and momentum equations read a consistent, immutable state to prepare time derivatives for the next iteration. In this normal case, it's OK to have the momentum equation after density.
However, now you are recalculating Vn in finally() from the ion velocity. This means that now finally() is inconsistent, with the density equation getting the "previous iteration".
If I'm right, this would still run but probably cause performance issues. I think you need to move the velocity calculation to the top.
| } | ||
|
|
||
| Field3D Nn_eq = (Nnlim * nu_cx + Ne * nu_rec) / softFloor(nu_cx + nu_iz, density_floor); | ||
| NVn = Nn_eq * U; |
There was a problem hiding this comment.
This is missing AA, so is off by a factor of 2 for d.
| } | ||
| } | ||
|
|
||
| Field3D Nn_eq = (Nnlim * nu_cx + Ne * nu_rec) / softFloor(nu_cx + nu_iz, density_floor); |
There was a problem hiding this comment.
There is a bug here, it floors collisionality with density_floor. The correct floor to use here is Rnn, which is derived from the "maximum neutral mean free path" neutral_lmax.
There was a problem hiding this comment.
I would also probably floor nu_cx and nu_rec to Rnn in the numerator for consistency, otherwise you'll have strange results.
| .doc("Allow anomalous diffusion into sheath?") | ||
| .withDefault<bool>(false); | ||
|
|
||
| anomalous_transport_momentum = options["anomalous_transport_momentum"] |
There was a problem hiding this comment.
What is this used for? The anomalous transport of momentum is just the viscosity eta. To disable it, just set eta to zero.
|
|
||
| } else if (passive_momentum) { | ||
| // NVn set from an equilibrium flow with the ions (not evolved): | ||
| // NVn = Nn_eq * Vi, Nn_eq = (Nn*nu_cx + Ne*nu_rec) / (nu_cx + nu_iz) |
There was a problem hiding this comment.
I haven't quite understood how this works, and why you aren't using plain Nn. Let's discuss this offline. It would be good to explain it in documentation. If these collisionalities are really needed, then I would ask that they are encapsulated in a function defined somewhere in the top of the file so they don't break up the flow so much. Maybe you could do this for the other collisionality block as well at that point. I can help with this if you need!
This is really my fault, I wrote this collisionality system and it's super ugly, so sorry... I will refactor it at some point!
| } | ||
| } | ||
|
|
||
| Field3D Nn_eq = (Nnlim * nu_cx + Ne * nu_rec) / softFloor(nu_cx + nu_iz, density_floor); |
There was a problem hiding this comment.
Equation is missing the pressure gradient term.
Also, could we rearrange this a bit for clarity? Calculating an intermediate Nn_eq makes it look like we are calculating density, which can be confusing on first read. I would prefer if we calculated NVn straight away, even if it's a long equation. Don't forget the mass factor as per the comment below.
| evolve_momentum = options["evolve_momentum"] | ||
| .doc("Evolve parallel neutral momentum?") | ||
| .withDefault<bool>(true); | ||
| passive_momentum = options["passive_momentum"] |
There was a problem hiding this comment.
I know that this name is now used also in the PR name, but I think there is an opportunity here to better communicate what this model is doing.
In effect, this model replaces the neutral momentum equation with a pressure-diffusion equation, just like what we use for cross-field transport: it's a balance of the neutral pressure gradient and momentum channels from the reactions (CX, IZ and REC).
If we called this option parallel_pressure_diffusion, then it would be clear that the pressure-diffusion model is used to replace the momentum equation. What do you think?
When I was reading passive_momentum I was always assuming this means Vn = Vi, because that's kind of what it sounds like to me.
… factor Three fixes to the quasi-static NVn calculation in passive_momentum mode: - Add missing parallel pressure-gradient term -(DnnNn/Pnlim)*Grad_par(Pn), the diffusive part of Bufferand 2024 eq (3) with D_n = Dnn/Tn. - Replace softFloor(nu_cx+nu_iz, density_floor) with nu_cx+nu_iz+Rnn: density_floor is dimensionally wrong as a collision-frequency floor; Rnn (neutral-neutral) is the physically motivated regularisation, consistent with how Dnn itself avoids division by zero. - Multiply both terms by AA: momentum sources in reaction.cxx use AA*velocity, so the quasi-static NVn = AA*[Nn_eq*U - D_n*Grad_par(Pn)]. Previously missing for deuterium (AA=2), harmless only for hydrogen (AA=1).
…ure_from guard Three changes from PR boutproject#592 review: - Rename parallel_pressure_diffusion (was passive_momentum): the option replaces the neutral momentum equation with a quasi-static pressure-diffusion balance, which better describes the physics than the old name. - evolve_momentum now defaults to !parallel_pressure_diffusion so users don't have to set it explicitly; raise BoutException if both are true simultaneously (was silently broken). - Raise BoutException if parallel_pressure_diffusion=true but temperature_from is not set: without an ion temperature source Pn stays zero, causing silent failures downstream.
…ure_from guard Three changes from PR boutproject#592 review: - Rename parallel_pressure_diffusion (was passive_momentum): the option replaces the neutral momentum equation with a quasi-static pressure-diffusion balance, which better describes the physics than the old name. - evolve_momentum now defaults to !parallel_pressure_diffusion so users don't have to set it explicitly; raise BoutException if both are true simultaneously (was silently broken). - Raise BoutException if parallel_pressure_diffusion=true but temperature_from is not set: without an ion temperature source Pn stays zero, causing silent failures downstream.
# Conflicts: # src/anomalous_diffusion.cxx # src/neutral_mixed.cxx
|
Update the comment according to the new pattern: PurposeAdd a reduced "parallel pressure-diffusion" neutral model to replace the earlier Change Summary
Validation
AI Assistance
Documentation
Review Notes
|
# Conflicts: # src/anomalous_diffusion.cxx
2 participants
Re-applies the previously-proposed neutral_passive_momentum feature onto current master (the original branch had drifted ~600 commits behind).
New options:
In passive_momentum mode NVn is set from an ion-equilibrium relation NVn = Nn_eq * Vi, with Nn_eq built from cx/iz/rec collision frequencies, and nu_cx/nu_iz/nu_rec/Nn_eq exposed as diagnostics.
The preconditioner gains a density-only branch (I - gammaDiv_perp(DnnGrad_perp))*ddt(Nn) when pressure is not evolved. Pressure/momentum terms, sources, diagnostics, time-derivative scaling, freeze_low_density and CHECK blocks are gated on the new flags. Reading another species' temperature declares the required GuardedOptions read permission.
Tested on a 3D MAST-U CD case: builds clean, runs to completion with finite fields, reproduces the original branch's passive behaviour (Tn tracks Ti to ~0.1%, 7 evolved fields, consistent neutral diagnostics).