Skip to content

Commit

Permalink
evolve_pressure: Important fix to P boundary
Browse files Browse the repository at this point in the history 
Must clear P parallel slices or Grad_par will use them.
Alternatively, setBoundaryTo should clear the parallel slices
if they are not set in the argument.
  • Loading branch information
@bendudson
bendudson committed Aug 17, 2024
1 parent dde009a commit 97a21fd
Show file tree
Hide file tree
Showing 2 changed files with 2 additions and 3 deletions.
4 changes: 1 addition & 3 deletions src/electron_force_balance.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,9 +14,7 @@ void ElectronForceBalance::transform(Options &state) {
throw BoutException("Cannot calculate potential and use electron force balance\n");
}

// Get the electron pressure
// Note: The pressure boundary can be set in sheath boundary condition
// which depends on the electron velocity being set here first.
// Get the electron pressure, with boundary condition applied
Options& electrons = state["species"]["e"];
Field3D Pe = GET_VALUE(Field3D, electrons["pressure"]);
Field3D Ne = GET_NOBOUNDARY(Field3D, electrons["density"]);
Expand Down
1 change: 1 addition & 0 deletions src/evolve_pressure.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -226,6 +226,7 @@ void EvolvePressure::finally(const Options& state) {

// Get updated pressure and temperature with boundary conditions
// Note: Retain pressures which fall below zero
P.clearParallelSlices();
P.setBoundaryTo(get<Field3D>(species["pressure"]));
Field3D Pfloor = floor(P, 0.0); // Restricted to never go below zero

Expand Down

3 comments on commit 97a21fd

@mikekryjak
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @bendudson, what does clearParallelSlices do?

@bendudson
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There are two ways to take parallel derivatives in BOUT++: One is to shift a field into field-aligned coordinates, and the other is to use other fields (the parallel slices) to store the values in the yup and ydown directions. When an operator like Grad_par operates on a field, it first checks to see if it has parallel slices. If it does then it uses them (no shifting needed), and if it doesn't then it does the to/from field aligned operation.

There are therefore two possible ways to apply boundary conditions:
(1) Shift field to aligned coordinates, apply the boundary, and shift back. The boundary values are stored in the guard cells of the field itself, not in its yup/down parallel slices. This is what Hermes-3 sheath_boundary etc. do.
(2) Store them in the parallel slices.

The sheath boundary components clear the parallel slices (so Grad_par etc won't use them), but in this case we keep a copy of the original pressure field in evolve_pressure, and set its boundary condition from the value in the state. This was done to help with handling of low-density regions. The sheath boundary condition cleared the parallel slices of the pressure in the state, but not for the copy stored in evolve_pressure. The result was that Grad_par used the wrong boundary condition on pressure.

I think the correct fix is probably to improve setBoundaryTo in BOUT++ (called here https://github.com/bendudson/hermes-3/blob/master/src/evolve_pressure.cxx#L219) so that it handles the presence (or lack) of poloidal slices.

@bendudson
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pull request with fix into BOUT++ next branch here: boutproject/BOUT-dev#2962

Please sign in to comment.