neutral_mixed MMS tests and expanded differential operator MMS tests
#327
Conversation
…or encountered: Field3D: Operation on empty data.
…e accumulation at isolated points on core boundary persists.
…. Strange accumulation at isolated points on core boundary persists." This reverts commit d7de1d3.
…r -- Error encountered: Field3D: Operation on empty data." This reverts commit 07a7b12.
…northog()." This reverts commit 690c1ff.
…pears to be 1 instead of 2.
… expected convergence order, permitting first-order accurate operators.
…t current test structure seems to require a different list and function call block in main.cxx due to differential operator functions have unique inputs, with a mix of const and non-const Field3D variables.
…g ddt for the model. This commit provides an initial job script and templates for the runs.
…llision frequency is probably not correctly captured.
…g easier. Put sensible default options into the test script.
… is small, increase length of timestep to give chance for larger error to develop.
…of these quantities.
…or in pressure diffusion.
This reverts commit 911ca20.
…t, make separate test directories for each test.
…tral_mixed. We test that the `ddt(Nd)`, `ddt(NVd)` and `ddt(Pd)` match the expected definitions from documentation, and we test that the RHS conserves `Nd` and `Pd` (or energy if `NVd` is evolved) by integrating `ddt(*)` over the volume. This test can be used to help give confidence when refactoring `neutral_mixed.cxx`. Note that the convergence order appears to be first-order accurate when `NVd` is evolved (or worse for energy conservation). This may be explained by the use of a slope limiter for terms involving `NVd`. Tests have been set to pass in the current state, but further investigations into the best numerical methods for preserving energy conservation should be carried out.
…at a momentum source can be specified in SI units for `neutral_mixed`.
…the hermes-3 jobs required.
|
The reason the CIs keep failing on this is that you're using the CVODE solver but SUNDIALS isn't available in your CI image. I've switched it over to PVODE, which seems to work. I've also updated the MMS testing script to print the output when there is an error, which should make such issues easier to diagnose in future. |
|
I'm struggling to work out why the CI is continuing to fail for Run 1947Error seems to occur after reading in and printing the second BOUT++ dataset for the final test case. Presumably it actually happens while reading the third (and final) dataset. Reported as a segfault. No backtrace is provided. Run 1948Error seemed to occurred when converting l2norm to a NumPy array for the final test case. This is after all of the datasets have been read in. Run 1949Same as Run 1947. Run 1950Error in final test case. I think this happens when turning the l2error into a NumPy array, again. Run 1955Same as 1947 and 1950, except this time there is a backtrace: Run 1956Again, seems to happen when converting the l2norm to a NumPy array for the final test case. This time there is no backtrace, just a low-level memory error. Some reflectionsThe evenly numbered runs are those labelled as being for the "pull request", while the odd number of runs are labelled as for the "push". In reality, both run exactly the same workflow. However, the "pull request" jobs fail when converting the l2norm to a NumPy array while the "push" jobs fail slightly earlier, while reading in the third dataset. The error messages that get printed vary from run to run, but not are particularly helpful in figuring out what is going wrong. |
|
Thanks @cmacmackin for your investigation so far. Maybe @ZedThree could chime in once he's back from leave - this seems particularly subtle.. |
| Dnn = (Tn / AA) / Rnn; | ||
| if (localstate.isSet("collision_frequency")) { | ||
| // Dnn = Vth^2 / sigma | ||
| Dnn = (Tn / AA) / (get<Field3D>(localstate["collision_frequency"]) + Rnn); |
There was a problem hiding this comment.
This is partly a question about the previous implementation, but when collision_frequency is set, we have collision_frequency + Rnn on the denominator: a) is that correct? b) should that be reflected in these changes?
We could do something like:
const Field3D Rnn = (rnn_override > 0.0) ? rnn_override : sqrt(floor(Tn, 1e-5) / AA) / neutral_lmax;
if (localstate.isSet("collision_frequency")) {
// Dnn = Vth^2 / sigma
Dnn = (Tn / AA) / (get<Field3D>(localstate["collision_frequency"]) + Rnn);
} else {
Dnn = (Tn / AA) / Rnn;
}
There was a problem hiding this comment.
Yes, this is correct. Rnn is a pseudo-collisionality which goes into the neutral diffusion model to limit its mean free path to a user set length, e.g. the machine size. It's a concept specific to this component and shouldn't live elsewhere.
There was a problem hiding this comment.
Should the new user-override also add on to collision_frequency?
There was a problem hiding this comment.
The user input is actually neutral_lmax, the max mean free path length scale. Rnn does not need to be a user input.
Both variables are poorly named, by the way. I have reworked this on my fork which improves the flux limitation for the neutral model which is due testing and merging once I get the time.
| print("../.././hermes_mms_tests -d "+workdir+" > "+workdir+"/output.txt") | ||
| os.system("../.././hermes_mms_tests -d "+workdir+" > "+workdir+"/output.txt") | ||
| s = os.system("../.././hermes_mms_tests -d "+workdir+" > "+workdir+"/output.txt") | ||
| if s != 0: | ||
| print(f"Command exited with status {s}. STDOUT printed below:") | ||
| with open(workdir + "/output.txt") as f: | ||
| print(f.read()) |
There was a problem hiding this comment.
Here's a slightly neater pattern:
| print("../.././hermes_mms_tests -d "+workdir+" > "+workdir+"/output.txt") | |
| os.system("../.././hermes_mms_tests -d "+workdir+" > "+workdir+"/output.txt") | |
| s = os.system("../.././hermes_mms_tests -d "+workdir+" > "+workdir+"/output.txt") | |
| if s != 0: | |
| print(f"Command exited with status {s}. STDOUT printed below:") | |
| with open(workdir + "/output.txt") as f: | |
| print(f.read()) | |
| workdir = pathlib.Path(workdir) | |
| command = f"../.././hermes_mms_tests -d {workdir}" | |
| print(command) | |
| s, out = launch(command, pipe=True) | |
| (workdir / "output.txt").write_text(out) | |
| if s != 0: | |
| print(f"Command exited with status {s}. STDOUT printed below:\n{out}") |
You can also use launch_safe to raise an exception if the command fails
(Also note to self: maybe launch/launch_safe should be able to write to a log file automatically?)
There was a problem hiding this comment.
This suggestion created error messages when I tried it locally. If you want to avoid using the os module, could you please make a PR with the appropriate changes?
P.S. errors were not due to failing to include pathlib, which I could manage, but due to something more mysterious
terminate called after throwing an instance of 'BoutExceptionterminate called after throwing an instance of 'BoutException'
what(): ERROR: Cannot split 60 Y points equally between 32 processors
| if not os.path.isdir(base_test_dir): | ||
| os.system("mkdir "+base_test_dir) | ||
| # create sub-directory, if required | ||
| if not sub_test_dir is None: | ||
| base_test_dir = base_test_dir + "/" + sub_test_dir | ||
| if not os.path.isdir(base_test_dir): | ||
| os.system("mkdir "+base_test_dir) |
There was a problem hiding this comment.
I know this is based on the other implementation, but I very much recommend using pathlib.Path for filesystem stuff in Python. This can a couple of lines:
| if not os.path.isdir(base_test_dir): | |
| os.system("mkdir "+base_test_dir) | |
| # create sub-directory, if required | |
| if not sub_test_dir is None: | |
| base_test_dir = base_test_dir + "/" + sub_test_dir | |
| if not os.path.isdir(base_test_dir): | |
| os.system("mkdir "+base_test_dir) | |
| base_test_dir = pathlib.Path(base_test_dir) | |
| base_test_dir.mkdir(parents=True, exist_ok=True) | |
| if sub_test_dir is not None: | |
| (base_test_dir / sub_test_dir).mkdir(parents=True, exist_ok=True) |
There was a problem hiding this comment.
Same as above, this suggestion created error messages when I tried it locally. If you want to avoid using the os module, could you please make a PR with the appropriate changes?
|
|
||
| return success, output_message | ||
|
|
||
| def run_neutral_mixed_manufactured_solutions_test(test_input): |
There was a problem hiding this comment.
What's actually different between this and run_manufactured_solutions_test? Is it just the base template input file? Maybe we could squash them together somehow?
There was a problem hiding this comment.
The run_manufactured_solutions_test() test calls a standalone executable to test the differential operators, run_neutral_mixed_manufactured_solutions_test() calls hermes-3 to test the neutral_mixed equations. I would not squash these functions, unless you can see a nice way to do it.
| const auto differential_operators_1_arg = { | ||
| nameandfunction1{"Div_par(f)", | ||
| [](const Field3D &f) { | ||
| return Div_par(f); | ||
| }}, | ||
| nameandfunction1{"Grad_par(f)", | ||
| [](const Field3D &f) { | ||
| return Grad_par(f); | ||
| }}, | ||
| }; |
There was a problem hiding this comment.
I think you should be able to do:
| const auto differential_operators_1_arg = { | |
| nameandfunction1{"Div_par(f)", | |
| [](const Field3D &f) { | |
| return Div_par(f); | |
| }}, | |
| nameandfunction1{"Grad_par(f)", | |
| [](const Field3D &f) { | |
| return Grad_par(f); | |
| }}, | |
| }; | |
| const auto differential_operators_1_arg = { | |
| nameandfunction1{"Div_par(f)", Div_par}, | |
| nameandfunction1{"Grad_par(f)", Grad_par}, | |
| }; |
There was a problem hiding this comment.
Attempting to do this leads to error messages like
error: cannot resolve overloaded function 'Div_par' based on conversion to type 'std::function<Field3D(const Field3D&)>'
57 | nameandfunction1{"Div_par(f)", Div_par},
In principle, I happy for this code to be simplified if possible, if you wanted to make a PR with working changes.
| for (const auto& difop: differential_operators_1_arg) { | ||
| if (difop.name.compare(differential_operator_name) == 0){ |
There was a problem hiding this comment.
I'm not 100% sure what's going on here. I think the intention is that we have a whole bunch of operators we want to iterate over but they take different arguments, which means they can't be stored in the same collection. So what we end up doing is having a bunch of different collections, and a loop for each one to both look for and call the particular operator.
If that's the case, then we can probably do something like use a base class with the common information, then we can have a single collection that we iterate over. We could also move the duplicated dump bits into a helper function. Something like:
class NamedFunctionBase {
std::string name;
std::string outname;
std::string expected_name;
public:
NamedFunctionBase(std::string name_, int i) :
name(std::move(name_)),
outname(fmt::format("result_{}", i),
expected_name(fmt::format("expected_result_{}", i)
{}
template<class... Args>
void operator()(Options& dump, Mesh& mesh, Args... args) {
dump[outname] = func(args...);
dump[outname].setAttributes({{"operator", name}});
Field3D expected_result{&mesh};
mesh.get(expected_result, expected_name, 0.0, false);
dump[expected_name] = expected_result;
}
};I think this is probably a little bit tricky, so I'm happy to have a chat about this, or have a stab at it myself.
There was a problem hiding this comment.
You have understood correctly -- initially there was a single collection of operators, but then different operator types were added, and this was my working solution. Happy for this to be simplified if it can be.
The changes included here aim to significantly expand the tests available in
tests/mms. The following features are added:tests/mms/perpendicular_laplacian.pymodule covering all of the differential operators that are present inneutral_mixed.nonorthogonal_test.pyandorthogonal_test.pyto test the newly added differential operators. This requires some changes tomain.cxx.neutral_mixed. We check that the definition of the right hand side matches the expected convergence order. Note that the expected convergence order is second order (2) whenevolve_momentum = false, but the observed convergence order is first order or worse whenevolve_momentum = true.neutral_mixed.neutral_mixed.The major addition is the test script
neutral_mixed_ddt_test.py, which runs the following four tests with prescribed functionsNd,Pd,NVdthat are constructed to be functions of(x,y)and satisfy the boundary conditions imposed inneutral_mixed.neutral_mixed/evolve_momentum_False_conservation_test_False. Check the definition of the RHS whenevolve_momentum = false, by calculating symbolicallyddt(Nd)andddt(Pd), and using these symbolic expressions as sources to cancel the numericalddt(Nd)andddt(Pd), such that the outputddt(Nd)andddt(Pd)are of order the discretisation error. The spatial resolution is scanned to check convergence, with the following results.neutral_mixed/evolve_momentum_False_conservation_test_True. Check that the RHS whenevolve_momentum = falseconservesNdandPd. The spatial resolution is scanned to check convergence, with the following results. Note that density is exactly conserved but pressure is only conserved within numerical error. A volume integration is done overddt(Nd)andddt(Pd)(without a symbolic source) to obtain these errors.neutral_mixed/evolve_momentum_True_conservation_test_False. Check the definition of the RHS whenevolve_momentum =true, by calculating symbolicallyddt(Nd),ddt(NVd)andddt(Pd), and using these symbolic expressions as sources to cancel the numericalddt(Nd),ddt(NVd)andddt(Pd), such that the outputddt(Nd),ddt(NVd)andddt(Pd)are of order the discretisation error. The spatial resolution is scanned to check convergence, with the following results.neutral_mixed/evolve_momentum_True_conservation_test_True. Check that the RHS whenevolve_momentum = trueconservesNdand energy. The spatial resolution is scanned to check convergence, with the following results. Note that conservation is first order accurate for density and less than first order for energy.These tests can be readily extended to test the
ddt()of other physics modules, or to 3D(x,y,z), allowing for rapid testing of the definition of the equations solved.Potential issues requiring further work:
neutral_mixeddocs and spurious (?) mass factor #328 for notes).neutral_mixedworking smoothly with nonorthogonal mesh in realistic physics case.