Start writing automated tests for MMS; renormalise manufactured distribution functions #64
Conversation
Relative errors mess up at points where the expected solution is 0.
Previously were always Float64.
Need loop that does not include vperp in the following code, so the previous begin_s_r_z_vperp_region() could cause errors.
Will allow it to be used in manufactured_solns.jl.
Makes manufactured dfni not exactly an even function in vpa, so that parallel velocity is not exactly zero and there is some z-advection.
Tricky to write down manufactured solution that exactly obeys "periodic" bc in vpa, so switch to "zero" bc instead. Add prefactor to manufactured f so it exactly obeys this bc at +/- Lvpa/2. Add odd component to manufactured f so upar will be non-zero (and non-constant). Done in a slightly strange way to ensure the additional component is odd, but never causes f to be negative anywhere. Tried adding non-zero upar by using a shifted Maxwellian type function, but this introduces erf() functions with non-constant arguments to the normalization correction (see below) which can be very expensive to evaluate, making the code slow. Ensure that the manufactured solution when integrated on -Lvpa/2<vpa<Lvpa/2 and 0<vperp<Lvperp is exactly equal to the manufactured n. This is necessary to ensure convergence, as the code only evaluates the numerical integral within these limits.
These points are not expected to agree with the manufactured rhs evaluation.
Normalise dfnn in manufactured solutions so that its integral on +/-Lv*/2 is exactly densn, for consistency with numerical integration routines. Add odd components to manufactured dfnn to avoid symmetry masking some bug. Rationalise argument order in manufactured solutions functions: generally r,z,vperp,vpa,vzeta,vr,vz except in the final float-evaluating function where the reverse vz,vr,vzeta,vpa,vperp,z,r is used to match the array index order.
Construct vrvzvzeta_dfni by substituting for vperp and vpa in terms of vzeta, vr, and vz. This means any normalisation, etc. that is applied to dfni is automatically carried over to vrvzvzeta_dfni. Similarly construct gav_dfnn by substituting for vzeta, vr, and vz in terms of vperp and vpa. Requires using a trick vzeta=vperp=-vr which effectively gyroaverages by getting rid of the particular form of odd components in vzeta and vr that were defined in dfnn_sym.
|
Regarding 1803b44 Thanks for the catch. I have also likely missed added sn (neutral species index) to parallelisation in several places, including this one. |
| ------- | ||
| (densi, phi, dfni) | ||
| """ | ||
| function manufactured_solutions_as_arrays( |
There was a problem hiding this comment.
Does a function like this make it possible to avoid passing the MMS functions around the main code?
If I not mistaken, there appears to be a speed cost associated with passing the runtimegenerated function -- it would be interesting to know if you notice a significant speedup by using this approach instead.
There was a problem hiding this comment.
I don't think there would be, mostly because this function allocates new arrays for its output, then returns them. Also it's re-building the manufactured-solutions functions each time it's called.
It should be possible to set up a struct in a way that everything in the main time-loop is fully typed. I'll have a look now and check if there is some type instability...
| # Because of particular form of odd component in dfnn_sym, if we set vzeta=-vr | ||
| # then we eliminate the odd component in the perpendicular velocities, | ||
| # effectively gyro-averaging | ||
| gav_dfnn = substitute(dfnn, Dict(vzeta=>vperp/sqrt(two), vr=>-vperp/sqrt(two), vz=>vpa)) |
There was a problem hiding this comment.
I am worried that this operation is only correct when the field line pitch is inifinite (or zero), so that z and b are codirectional. I think that you are wanting vzeta and vr to be the velocities in the plane perpendicular to the magnetic field line, which they are not. In my original 'gyroaveraged' functions code, I was exploiting that vperp^2 + vpa^2 = vtot^2 = vz^2 + vr^2 + vzeta^2, so the gyroaverage was trivial. Am I missing the point here?
There was a problem hiding this comment.
I made a mistake here - was thinking z and b are in the same direction. Was trying to be clever and give non-zero flow, but this doesn't work 😞
|
Looking at the manufactured_solns.jl file, I have two immediate comments. You say in the commit message that the need for the boundary terms goes away as soon as the length of the domains are large enough. Since computing the boundary terms is painful and not automated, is it possible to dispense with them by making sure that the tests cover a reasonably large range in velocities? I am not sure that the method that you are using for the boundary terms will carry over to the wall boundary test currently in the vanilla branch. For the neutrals, the distribution function must be proportional to the Knudsen distribution in order to match the boundary condition. Finally, a matter of taste. I rather liked having the kinetic equation and the drift-kinetic equation written on one line, as we would write it for a slide. I thought that this would have value by the way of self documentation. I can see why it could be convenient to add the terms one by one to a "rhs". If we want to go that way I think we need to make sure that there are adequate comments documenting which equation is being solved in the source code near this function. Since Julia permits expressions extended over multiple lines, I would have thought that we could avoid spreading the equation over too much of the source code. |
Just type-parameterise the struct, replacing the NCDatasets-dependent variable types which can be broken by changes in NCDatasets.jl version, etc.
johnomotani
force-pushed
the
radial-vperp-standard-DKE-auto-MMS
branch
from
b3349fe
to
2fae0b8
Compare
...as well as cases with neutrals. Without neutrals the code has lower dimensionality, so it is practical to test at higher resolution (before running out of memory).
Only increase n_element by 2x for velocity dimensions compared to spatial dimensions, instead of 4x. Makes tests cheaper.
Start on automated tests using MMS. Time-dependent test
test/manufactured_solutions_tests.jlis in very rough state. Another test which only tests the 'rhs' evaluation at t=0,test/manufactured_ddt_tests.jlis closer to being finished, but still just prints the errors, convergence orders (and a bit of other debugging info) - still needs a decision on what resolution parameters to test automatically, and a pass/fail check and tolerance setting. Note that increasing the resolution too much can easily lead to running out of RAM due to the 5d neutral distribution function.To get the MMS tests to converge nicely, redefine the manufactured distribution functions so that the ion distribution function when integrated over 0<=vperp<=Lvperp, -Lvpa/2<=vpa<=Lvpa/2 gives exactly the manufactured ion density, and the neutral distribution function integrated on -Lvzeta/2<=vzeta<=Lvzeta/2, -Lvr/2<=vr<=Lvr/2, -Lvz/2<=vz<=Lvz/2 gives exactly the manufactured neutral density. This avoids the error tending to a floor related to the Lv* settings (although this error can be made very small - below machine precision - by increasing Lv* by a factor of 2 or so).
The charge exchange term still only converges at somewhere between 1st and 2nd order (for the resolutions I have tested) - I think this is expected due to the use of linear interpolation in this term. This would probably be improved dramatically by @mrhardman's suggestion of using a vpa,vperp,gyrophase grid for the neutrals.
Fixes a couple of small bugs: copy-paste error in an index in charge_exchange_3V!() function, could have had a significant impact (db2cf67); wrong (?) region in
ssp_rk!()may cause errors for some combinations of dimension parallelisation 1803b44.