diamagnetic_drift: Boundary flux changes - #539
Conversation
Previous method, using a varying weighting between divergence and gradient form, may not have preserved fluxes. Introducing a new method that averages over the core boundary and prevents fluxes through boundaries by default.
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #539 +/- ##
==========================================
+ Coverage 49.09% 49.44% +0.34%
==========================================
Files 96 97 +1
Lines 10037 10102 +65
Branches 1452 1466 +14
==========================================
+ Hits 4928 4995 +67
+ Misses 4604 4591 -13
- Partials 505 516 +11 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
Split the transform function into pieces and add unit tests. Includes checks that the operator conserves fluxes in both diamagnetic and gradient form.
bendudson
changed the title
No longer blending divergence and gradient forms.
cmacmackin
left a comment
cmacmackin
left a comment
There was a problem hiding this comment.
I don't like the way public functions have been introduced in this. They are not necessary to test the behaviour of the component and will increase code maintenance burdens.
Even if you decide to keep calculateDivergenceForm and calculateGradientForm public, I would strongly urge you to make addDiamagneticDrift private. I see no value in calling it directly in your unit tests and doing so requires setting up GuardedOptions objects. Doing that will require you to keep track of the permissions set by DiamagneticDrift and make your tests significantly more fragile as a result..
|
|
||
| divergence_form = options["divergence_form"] | ||
| .doc("Use divergence form of diamagnetic drift?") | ||
| .withDefault<bool>(true); |
There was a problem hiding this comment.
Call this option "form" and use BoutEnum to provide options with meaningful names (e.g., divergence and gradient).
| if (!average_core) { | ||
| throw BoutException("DiamagneticDrift::coreAverage requires average_core=true"); | ||
| } | ||
|
|
There was a problem hiding this comment.
We can avoid this check if we make coreAverage is a standalone function that take core_ring_volume and cell_volume as arguments.
| if (fabs(charge) < 1e-5) { | ||
| throw BoutException( | ||
| "DiamagneticDrift::calculateGradientForm requires a non-zero charge"); | ||
| } |
There was a problem hiding this comment.
Don't need this, as we already check the charge before calling.
| if (fabs(charge) < 1e-5) { | ||
| throw BoutException( | ||
| "DiamagneticDrift::calculateDivergenceForm requires a non-zero charge"); | ||
| } |
There was a problem hiding this comment.
Don't need this, as we already check the charge before calling.
| for (int jx = mesh->xstart; jx <= mesh->xend; ++jx) { | ||
| for (int jy = mesh->ystart; jy <= mesh->yend; ++jy) { | ||
| for (int jz = mesh->zstart; jz <= mesh->zend; ++jz) { | ||
| result += cell_volume(jx, jy) * field(jx, jy, jz); | ||
| } | ||
| } | ||
| } |
There was a problem hiding this comment.
| for (int jx = mesh->xstart; jx <= mesh->xend; ++jx) { | |
| for (int jy = mesh->ystart; jy <= mesh->yend; ++jy) { | |
| for (int jz = mesh->zstart; jz <= mesh->zend; ++jz) { | |
| result += cell_volume(jx, jy) * field(jx, jy, jz); | |
| } | |
| } | |
| } | |
| BOUT_FOR(i, mesh->getRegion3D("RGN_NOBNDRY")) { | |
| result += cell_volume[i] * field[i]; | |
| } |
| FakeMesh div_mesh(6, 5, 7, *bout::globals::mpi); | ||
| div_mesh.xstart = 2; | ||
| div_mesh.xend = div_mesh.LocalNx - 3; | ||
| div_mesh.createDefaultRegions(); | ||
| div_mesh.setCoordinates(nullptr); | ||
| auto div_coords = makeUnitCoordinates(&div_mesh); | ||
| div_mesh.setCoordinates(div_coords); | ||
| div_mesh.setGridDataSource(new FakeGridDataSource()); | ||
| div_mesh.createBoundaryRegions(); | ||
|
|
||
| ScopedGlobalMesh use_div_mesh(&div_mesh); |
There was a problem hiding this comment.
As mentioned above, consider moving this logic into a fixture, to keep the body of the test cleaner.
| bndry_flux = options["bndry_flux"] | ||
| .doc("Allow fluxes through boundary?") | ||
| .withDefault<bool>(false); |
There was a problem hiding this comment.
Seeing as we're making breaking changes anyway, we might as well rename this "boundary_flux".
| }); | ||
| Field3D temperature{1.0}; | ||
|
|
||
| Field3D sink = component.calculateDivergenceForm(quantity, temperature, 1.0); |
There was a problem hiding this comment.
As described previously, I think it is better to do this with calls to DiamagneticDirft::transform. If you really want to do it this way, however, there is no need to set average_core or diveregence_form.
| ASSERT_NEAR(volumeIntegral(sink), 0.0, 1e-12); | ||
| } | ||
|
|
||
| TEST_F(WideDiamagneticDriftTest, GradientFormIsConservativeWithConstantBoundaryProduct) { |
There was a problem hiding this comment.
What is meant by BoundaryProduct in this test name. I see the constant boundary, but where is the product?
| a flux through the boundary, while the gradient form can be used to avoid | ||
| those unphysical boundary fluxes. | ||
|
|
||
| In Hermes-3 these are selected with the boolean option ``divergence_form``: |
There was a problem hiding this comment.
Use an enum for this, instead of a boolean.
2 participants
Purpose
Ensure that the divergence of the diamagnetic drift is zero, and that fluxes through the core boundary doesn't result in unphysical boundary layers.
Change Summary
transform_implfunction into public functionsValidation
Unit tests of steps in the calculation, following refactor. Includes tests of the conservation properties.
AI Assistance
Codex/GPT-5.4 wrote most of the code. I planned, reviewed and adjusted the changes.
Documentation
The diamagnetic drift equation documentation has been updated.
Review Notes