Skip to content

Commit

Permalink
Fixing bug in calculation of non-uniform correction
Browse files Browse the repository at this point in the history
As pointed out by Haruki Seto, the non-uniform correction
terms d1_dx and d1_dy were not recalculated if dx and dy were
changed.

This fix moves the calculation of d1_dx and d1_dy into the
geometry() function, which should be called whenever metric
tensor components are changed.
  • Loading branch information
bendudson committed Nov 6, 2017
1 parent a58b9ab commit ddfc374
Showing 1 changed file with 21 additions and 19 deletions.
40 changes: 21 additions & 19 deletions src/mesh/coordinates.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -169,26 +169,7 @@ Coordinates::Coordinates(Mesh *mesh) {
throw BoutException("Differential geometry failed\n");
}

//////////////////////////////////////////////////////
/// Non-uniform meshes. Need to use DDX, DDY

OPTION(Options::getRoot(), non_uniform, false);

Field2D d2x, d2y; // d^2 x / d i^2
// Read correction for non-uniform meshes
if (mesh->get(d2x, "d2x")) {
output_warn.write("\tWARNING: differencing quantity 'd2x' not found. Calculating from dx\n");
d1_dx = mesh->indexDDX(1. / dx); // d/di(1/dx)
} else {
d1_dx = -d2x / (dx * dx);
}

if (mesh->get(d2y, "d2y")) {
output_warn.write("\tWARNING: differencing quantity 'd2y' not found. Calculating from dy\n");
d1_dy = mesh->indexDDY(1. / dy); // d/di(1/dy)
} else {
d1_dy = -d2y / (dy * dy);
}

if (mesh->get(ShiftTorsion, "ShiftTorsion")) {
output_warn.write("\tWARNING: No Torsion specified for zShift. Derivatives may not be correct\n");
Expand Down Expand Up @@ -365,6 +346,27 @@ int Coordinates::geometry() {

mesh->communicate(com);

//////////////////////////////////////////////////////
/// Non-uniform meshes. Need to use DDX, DDY

OPTION(Options::getRoot(), non_uniform, false);

Field2D d2x, d2y; // d^2 x / d i^2
// Read correction for non-uniform meshes
if (mesh->get(d2x, "d2x")) {
output_warn.write("\tWARNING: differencing quantity 'd2x' not found. Calculating from dx\n");
d1_dx = mesh->indexDDX(1. / dx); // d/di(1/dx)
} else {
d1_dx = -d2x / (dx * dx);
}

if (mesh->get(d2y, "d2y")) {
output_warn.write("\tWARNING: differencing quantity 'd2y' not found. Calculating from dy\n");
d1_dy = mesh->indexDDY(1. / dy); // d/di(1/dy)
} else {
d1_dy = -d2y / (dy * dy);
}

return 0;
}

Expand Down

0 comments on commit ddfc374

Please sign in to comment.