From ddfc3742396062609ec3d30a983270db9882a649 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Mon, 6 Nov 2017 10:42:34 +0000 Subject: [PATCH] Fixing bug in calculation of non-uniform correction 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. --- src/mesh/coordinates.cxx | 40 +++++++++++++++++++++------------------- 1 file changed, 21 insertions(+), 19 deletions(-) diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 485ea4de5a..eb2ef45eb8 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -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"); @@ -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; }