Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Simplify calculation of coefficients in moment-kinetic equation #320

Draft
wants to merge 15 commits into
base: master
Choose a base branch
from

Conversation

johnomotani
Copy link
Collaborator

@johnomotani johnomotani commented Feb 13, 2025

Saving the time derivatives of the moments as they are calculated for the time-advance of the moments, and then using those time derivatives for the coefficients in the moment-kinetic equation, simplifies the code, making fewer places that need to be updated if new
terms are added (e.g. if new terms are added to the RHS of the drift-kinetic equation, only the moment equations and the RHS of the moment-kinetic equation change, not the $\dot z$, $\dot w_\parallel$, $\dot F$ coefficients).

Updates expected values in some tests: discretization errors are changed slightly by the rearrangement of the moment kinetic advection coefficients, and the tests are generally poorly resolved (in order to be fast), so are sensitive to this.

There is a question whether making this change is a good idea, with the concern that there may be some desirable cancellations between terms that would take place analytically but might not happen exactly numerically. Possibly the most concerning is the cancellation of the parallel electric field, but I think that one should be exact up to machine precision, as

dupar_dt = dnupar_dt / 	n + ...
                 = (0.5 * Ez * n + ...) / n + ...
wpa_speed = (0.5 * Ez - dupar_dt + ...) / vth

So that we end up evaluating numerically 0.5 * Ez - 0.5 * Ez * n / n - there will be some rounding error, but no discretization error (for other terms there might be a discretization error due to something like d(n*upar)/dz - n*dupar/dz - upar*dn/dz). For the corresponding lines of code, see:

dupar_dt[iz,ir,is] = (dnupar_dt[iz,ir,is] - upar[iz,ir,is] * dn_dt[iz,ir,is]) / n[iz,ir,is]

0.5*geometry.bzed[iz,ir]*fields.Ez[iz,ir]*density[iz,ir,is])

(0.5 * Ez[iz,ir]
- (dupar_dt[iz,ir,is] + (vth[iz,ir,is] * wpa + upar[iz,ir,is]) * dupar_dz[iz,ir,is])
- wpa * (dvth_dt[iz,ir,is] + (vth[iz,ir,is] * wpa + upar[iz,ir,is]) * dvth_dz[iz,ir,is])
) / vth[iz,ir,is]

This will allow the time derivatives of moments to be used directly in
the kinetic equation.
Simplifies code, making fewer places that need to be updated if new
terms are added (now only RHS of kinetic equation and moment equations,
not zdot, wpadot, Fdot coefficients).

Updates expected values in some tests: discretization errors are changed
by the rearrangement of the moment kinetic advection coefficients, and
the tests are generally poorly resolved (in order to be fast), so are
sensitive to this.
@johnomotani johnomotani requested a review from mabarnes February 13, 2025 17:36
@mrhardman
Copy link
Collaborator

mrhardman commented Feb 14, 2025

I am in favour of anything which simplifies the upgrading of the model. Thanks for looking at the possibility that this might work!

Looking in the code it seems that you have $\partial n / \partial t$, $\partial u / \partial t$, and $\partial p / \partial t$ left explicitly in place everywhere, so the equations are now written as the would be if you normalise them in the moment kinetic way but don't substitute for $\partial n / \partial t$, $\partial u / \partial t$, and $\partial p / \partial t$ -- Is this what you mean when you say we no longer have to change $\dot{z}$, $\dot{w_\|}$, and $\dot{F}$ when we update the DKE and the moment eqns?

Is there a way to make things even more general, so that the function which computes $\dot{z}$, $\dot{w_\|}$, and $\dot{F}$ for the case where no moments are split can be used as a starting point to add the moment kinetic corrections to? This way, we would not have to implement the $E \times B$ and magnetic drifts, e.g., in multiple places.

EDIT: I suppose here the challenge is that $w_\|$, $w_\perp$ have different meanings in the moment-split and standard implementations -- but we could probably write a function to compute the drifts that took as input the "true" velocity $v_\|$ and $v_\perp$, which can easily be computed from the normalised values. These functions can be called in each of the $\dot{z}$, $\dot{w_\|}$, and $\dot{F}$ implementations, so at least the drift implementation is in a single place.

We don't currently use a time advance scheme where the electron density
(equal to the ion density) is updated during the same step as the
electron pressure (because explicit and implicit steps are separate even
if the electrons are updated directly on the ion timestep), but if we
did, and also used the temperature equation option for the electrons
(instead of the default pressure equation), then the implementation
would be wrong because it did not distinguish the density at the
beginning and the end of the timestep (compare the implementation of the
parallel velocity update using a momentum equation in force_balance.jl).
This commit fixes this potential error, by passing in separate
electron_density_out and electron_density in to
electron_energy_equation!().
@johnomotani
Copy link
Collaborator Author

@mrhardman ah, my comment was not correct (I'll edit to try and clarify). I was only thinking of adding new terms to the RHS of the kinetic equation, not potentially adding drift terms or other modifications to the free-streaming or acceleration. Potentially we could do something along the lines you're suggesting for drift terms - I think it's something to come back to when we're implementing the 2D version of moment-kinetics.

Yes, using arrays with dn/dt, etc. instead of substituting by hand using the moment equations, then evaluating the resulting expressions is what this PR does.

@johnomotani johnomotani force-pushed the simplify-moment-kinetic-advection-coefficients branch from 8d9688a to 4afa758 Compare February 16, 2025 13:46
@mrhardman
Copy link
Collaborator

@mrhardman ah, my comment was not correct (I'll edit to try and clarify). I was only thinking of adding new terms to the RHS of the kinetic equation, not potentially adding drift terms or other modifications to the free-streaming or acceleration. Potentially we could do something along the lines you're suggesting for drift terms - I think it's something to come back to when we're implementing the 2D version of moment-kinetics.

Yes, using arrays with dn/dt, etc. instead of substituting by hand using the moment equations, then evaluating the resulting expressions is what this PR does.

Yes, in the best case scenario when modifying the advection velocities and accelerations we would have to modify both $\dot{z}$, $\dot{w}_\|$, $\dot{w}_\perp$, and $\partial n/ \partial t$, $\partial u_\|/ \partial t$, and $\partial p/ \partial t$. I think my suggestion was aimed at keeping the effort in those modifications to the theoretical minimum, happy for you to consider this at a later date if that is not the point of this PR.

LucasMontoya4 and others added 6 commits February 17, 2025 14:26
…n-coefficients' into simplify-moment-kinetic-advection-coefficients
…imes, rather than recursively updating each element in the time derivative calculator in get_variable().
By fixing _get_fake_moments_fields_scratch(), we can apply the fix in
a single place, instead of everwhere _get_fake_moments_fields_scratch()
is called.
…etic-advection-coefficients' into simplify-moment-kinetic-advection-coefficients
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants