Skip to content

Preview: Use cell area and volume for finite volume code#3392

Open
dschwoerer wants to merge 35 commits into
move-opsfrom
move-ops-area
Open

Preview: Use cell area and volume for finite volume code#3392
dschwoerer wants to merge 35 commits into
move-opsfrom
move-ops-area

Conversation

@dschwoerer

Copy link
Copy Markdown
Contributor

A pr of #3326 into #3335 to show the facility of cell volumes and areas

@dschwoerer

Copy link
Copy Markdown
Contributor Author

@github-actions github-actions Bot left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

clang-tidy made some suggestions

Comment thread include/bout/fv_ops.hxx
Comment thread include/bout/fv_ops.hxx Outdated
Comment thread src/mesh/coordinates.cxx
ASSERT0(mesh->xstart > 0);
BOUT_FOR(i, area_centre.getRegion("RGN_NOX")) {
(*_cell_area_xlow)[i] = 0.5 * (area_centre[i] + area_centre[i.xm()]);
(*_cell_area_xhigh)[i] = 0.5 * (area_centre[i] + area_centre[i.xp()]);

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

The factor 0.5 does imply here that the position of the interface is right at the geometric middle of the cells right? I don't see a reason to not weight this by the actual size of the cell.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Right now it moves over what the FA code does. I think hypnotoad might have more details about the size and direction of the actual boundary of the cell, but unless hypnotoad adds something like the dagp coefficients, like zoidberg does, it will not be possible to correctly reconstruct it.

So 0.5 is motivated, by several facts

  • it is easy
  • it is in line with the current FA code
  • Anything more complicated immediately needs to go to quite a level of detail, to include extending in the x direction, but also includes shearing, so it needs to depends on g_11, g_13 and g_12

If you have a formulation that is numerically stable, and can show that it gives better results, I am happy to change it. But to me that feels like overkill.

Comment thread src/mesh/coordinates.cxx
}
BOUT_FOR(i, By_c.getRegion("RGN_NOY")) {
jxz_ylow[i] = By_c[i] / By_l[i] * jxz_centre[i];
jxz_yhigh[i] = By_c[i] / By_h[i] * jxz_centre[i];

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

@dschwoerer We talked about this, again I am not sure that we are taking into account the "proper" conservation of magnetic flux, By here is just the toroidal field and not the total field strength. Maybe @bendudson or @bshanahan have an idea on how to handle this properly. I would argue that when the field line travels into an area of the same toroidal but a different poloidal field, this effect should have an impact on the cell area but here it does not. I know that this is tricky to do because right now we assume that perp = poloidal.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Lets iterate the same arguments another round :-)

By is the part that is orthogonal to the plane.
If we want to follow the magnetic flux tube, we need to take the area orthogonal to the magnetic field.
That area in the RZ cross section area_{rz} - and the area orthogonal to the magnetic field is smaller by a factor of By / Bmag, a sign of By changes the direction of the area. Thus the total flux would be Bmag * By / Bmag * area_xz = By * area_xz.

I do not see any further assumptions needed for using this expressions.

Certainly, we make further assumptions, but I think if they need to be corrected here in this expressions, we should just remove those assumptions, rather then making it wrong here, to correct for some assumptions some where else.

Note that this is also correct within the assumption that By = Bmag.

Comment thread src/mesh/coordinates.cxx
BOUT_FOR(i, area_centre.getRegion("RGN_NOZ")) {
(*_cell_area_zlow)[i] = 0.5 * (area_centre[i] + area_centre[i.zm()]);
(*_cell_area_zhigh)[i] = 0.5 * (area_centre[i] + area_centre[i.zp()]);
}

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Same argument as before

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Same as above.

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