From e9e00088f635e4c33e05e8172f7d21e33d0caeff Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 14 May 2018 21:40:25 +0100 Subject: [PATCH 1/4] Remove mixed field-aligned/non-field-aligned derivatives Advective and flux derivatives, which take v and f as inputs, had cases using yup/ydown fields for one of v and f, but not the other. These do not make sense as multiplication of a field in field-aligned coordinates with another in non-field-aligned coordinates is incorrect. Delete last of these cases and fall back to converting both v and f to field-aligned if either does not have yup/ydown fields. Also add some comments where mixed cases were removed before. --- src/mesh/difops.cxx | 57 ++++----------------------------------- src/mesh/index_derivs.cxx | 22 +++++++-------- 2 files changed, 15 insertions(+), 64 deletions(-) diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index 58625c77d6..0199d991a7 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -342,60 +342,13 @@ const Field3D Vpar_Grad_par_LCtoC(const Field3D &v, const Field3D &f, REGION reg } } } - else if (vUseUpDown) { - // Only v has up/down fields - // f must shift to field aligned coordinates - Field3D f_fa = vMesh->toFieldAligned(f); - - BOUT_OMP(parallel) { - stencil fval, vval; - BOUT_FOR_INNER(i, vMesh->getRegion3D(region_str)) { - fval.mm = f_fa[i.ymm()]; - fval.m = f_fa[i.ym()]; - fval.c = f_fa[i]; - fval.p = f_fa[i.yp()]; - fval.pp = f_fa[i.ypp()]; - - vval.m = v.ydown()[i.ym()]; - vval.c = v[i]; - vval.p = v.yup()[i.yp()]; - - // Left side - result[i] = (vval.c >= 0.0) ? vval.c * fval.m : vval.c * fval.c; - // Right side - result[i] -= (vval.p >= 0.0) ? vval.p * fval.c : vval.p * fval.p; - } - } - } - else if (fUseUpDown) { - // Only f has up/down fields - // v must shift to field aligned coordinates - Field3D v_fa = vMesh->toFieldAligned(v); - - BOUT_OMP(parallel) { - stencil fval, vval; - BOUT_FOR_INNER(i, vMesh->getRegion3D(region_str)) { - fval.m = f.ydown()[i.ym()]; - fval.c = f[i]; - fval.p = f.yup()[i.yp()]; - - vval.mm = v_fa[i.ymm()]; - vval.m = v_fa[i.ym()]; - vval.c = v_fa[i]; - vval.p = v_fa[i.yp()]; - vval.pp = v_fa[i.ypp()]; - - // Left side - result[i] = (vval.c >= 0.0) ? vval.c * fval.m : vval.c * fval.c; - // Right side - result[i] -= (vval.p >= 0.0) ? vval.p * fval.c : vval.p * fval.p; - } - } - } else { // Both must shift to field aligned - Field3D v_fa = vMesh->toFieldAligned(v); - Field3D f_fa = vMesh->toFieldAligned(f); + // (even if one of v and f has yup/ydown fields, it doesn't make sense to + // multiply them with one in field-aligned and one in non-field-aligned + // coordinates) + Field3D v_fa = mesh->toFieldAligned(v); + Field3D f_fa = mesh->toFieldAligned(f); BOUT_OMP(parallel) { stencil fval, vval; diff --git a/src/mesh/index_derivs.cxx b/src/mesh/index_derivs.cxx index d872b70ab4..56db9871eb 100644 --- a/src/mesh/index_derivs.cxx +++ b/src/mesh/index_derivs.cxx @@ -2158,12 +2158,8 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo func = lookupFunc(table, method); } - // There are four cases, corresponding to whether or not f and v - // have yup, ydown fields. - - // If vUseUpDown is true, field "v" has distinct yup and ydown fields which - // will be used to calculate a derivative along - // the magnetic field + // If *UseUpDown is true, field "*" has distinct yup and ydown fields which + // will be used to calculate a derivative along the magnetic field bool vUseUpDown = (v.hasYupYdown() && ((&v.yup() != &v) || (&v.ydown() != &v))); bool fUseUpDown = (f.hasYupYdown() && ((&f.yup() != &f) || (&f.ydown() != &f))); @@ -2195,6 +2191,9 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo } } else { // Both must shift to field aligned + // (even if one of v and f has yup/ydown fields, it doesn't make sense to + // multiply them with one in field-aligned and one in non-field-aligned + // coordinates) Field3D v_fa = this->toFieldAligned(v); Field3D f_fa = this->toFieldAligned(f); BOUT_OMP(parallel) { @@ -2788,12 +2787,8 @@ const Field3D Mesh::indexFDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo Field3D result(this); result.allocate(); // Make sure data allocated - // There are four cases, corresponding to whether or not f and v - // have yup, ydown fields. - - // If vUseUpDown is true, field "v" has distinct yup and ydown fields which - // will be used to calculate a derivative along - // the magnetic field + // If *UseUpDown is true, field "*" has distinct yup and ydown fields which + // will be used to calculate a derivative along the magnetic field bool vUseUpDown = (v.hasYupYdown() && ((&v.yup() != &v) || (&v.ydown() != &v))); bool fUseUpDown = (f.hasYupYdown() && ((&f.yup() != &f) || (&f.ydown() != &f))); @@ -2830,6 +2825,9 @@ const Field3D Mesh::indexFDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo } } else { // Both must shift to field aligned + // (even if one of v and f has yup/ydown fields, it doesn't make sense to + // multiply them with one in field-aligned and one in non-field-aligned + // coordinates) Field3D v_fa = this->toFieldAligned(v); Field3D f_fa = this->toFieldAligned(f); BOUT_OMP(parallel) { From 90648ba8a4d7cafb248807f81300df8c9fa3e081 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 14 May 2018 21:31:39 +0100 Subject: [PATCH 2/4] Add missing fromFieldAligned calls In several y-derivatives where we have to convert the input to field-aligned coordinates, the result was not transformed back to the original coordinates. Also remove uses of global 'mesh' in Vpar_Grad_par_LCtoC() --- src/mesh/difops.cxx | 6 ++++-- src/mesh/index_derivs.cxx | 4 ++++ 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index 0199d991a7..8b8ff649ce 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -347,8 +347,8 @@ const Field3D Vpar_Grad_par_LCtoC(const Field3D &v, const Field3D &f, REGION reg // (even if one of v and f has yup/ydown fields, it doesn't make sense to // multiply them with one in field-aligned and one in non-field-aligned // coordinates) - Field3D v_fa = mesh->toFieldAligned(v); - Field3D f_fa = mesh->toFieldAligned(f); + Field3D v_fa = vMesh->toFieldAligned(v); + Field3D f_fa = vMesh->toFieldAligned(f); BOUT_OMP(parallel) { stencil fval, vval; @@ -366,6 +366,8 @@ const Field3D Vpar_Grad_par_LCtoC(const Field3D &v, const Field3D &f, REGION reg // Right side result[i] -= (vval.p >= 0.0) ? vval.p * fval.c : vval.p * fval.p; } + + result = vMesh->fromFieldAligned(result); } } diff --git a/src/mesh/index_derivs.cxx b/src/mesh/index_derivs.cxx index 56db9871eb..886788ab3d 100644 --- a/src/mesh/index_derivs.cxx +++ b/src/mesh/index_derivs.cxx @@ -2224,6 +2224,8 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo result[i] = func(vval, fval); } } + + result = this->fromFieldAligned(result); } } else { // Non-staggered case @@ -2860,6 +2862,8 @@ const Field3D Mesh::indexFDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo result[i] = func(vval, fval); } } + + result = this->fromFieldAligned(result); } result.setLocation(outloc); From 54494d4cfcfb99f17b66963df80967da9b2285a2 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 3 Jul 2018 09:58:48 +0100 Subject: [PATCH 3/4] Fix f->f_fa Was previously a typo that resulted in using non-field-aligned f in VDDY, which is incorrect. --- src/mesh/index_derivs.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mesh/index_derivs.cxx b/src/mesh/index_derivs.cxx index 886788ab3d..c05991520d 100644 --- a/src/mesh/index_derivs.cxx +++ b/src/mesh/index_derivs.cxx @@ -2207,7 +2207,7 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo fval.mm = f_fa[i.ymm()]; fval.m = f_fa[i.ym()]; - fval.c = f[i]; + fval.c = f_fa[i]; fval.p = f_fa[i.yp()]; fval.pp = f_fa[i.ypp()]; From 1342d874cbc3b00e6fd6fb64c1899d1ecb436a5c Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 3 Jul 2018 12:16:41 +0100 Subject: [PATCH 4/4] Set location of result as soon as it is declared in index_derivs.cxx The location of the result will in future need to be set before shifting it from field-aligned coordinates, if this is necessary. So it is safer to set the location to outloc as soon as 'Field3D result' is declared. Also in a couple of places just 'return apply*diff(...)' instead of creating an unneeded intermediate variable 'result'. --- src/mesh/index_derivs.cxx | 70 +++++++++++---------------------------- 1 file changed, 20 insertions(+), 50 deletions(-) diff --git a/src/mesh/index_derivs.cxx b/src/mesh/index_derivs.cxx index c05991520d..f684cce292 100644 --- a/src/mesh/index_derivs.cxx +++ b/src/mesh/index_derivs.cxx @@ -662,6 +662,7 @@ const Field2D Mesh::applyXdiff(const Field2D &var, Mesh::deriv_func func, Field2D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); if (this->StaggerGrids && (outloc != inloc)) { // Staggered differencing @@ -752,8 +753,6 @@ const Field2D Mesh::applyXdiff(const Field2D &var, Mesh::deriv_func func, } } - result.setLocation(outloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false; @@ -787,6 +786,7 @@ const Field3D Mesh::applyXdiff(const Field3D &var, Mesh::deriv_func func, Field3D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); if (this->StaggerGrids && (outloc != inloc)) { // Staggered differencing @@ -877,8 +877,6 @@ const Field3D Mesh::applyXdiff(const Field3D &var, Mesh::deriv_func func, } } - result.setLocation(outloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false; @@ -912,6 +910,7 @@ const Field2D Mesh::applyYdiff(const Field2D &var, Mesh::deriv_func func, CELL_L Field2D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); if (this->ystart > 1) { // More than one guard cell, so set pp and mm values @@ -943,8 +942,6 @@ const Field2D Mesh::applyYdiff(const Field2D &var, Mesh::deriv_func func, CELL_L } } - result.setLocation(outloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_yup = result.bndry_ydown = false; @@ -977,6 +974,7 @@ const Field3D Mesh::applyYdiff(const Field3D &var, Mesh::deriv_func func, CELL_L Field3D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); if (var.hasYupYdown() && ((&var.yup() != &var) || (&var.ydown() != &var))) { // Field "var" has distinct yup and ydown fields which @@ -1122,8 +1120,6 @@ const Field3D Mesh::applyYdiff(const Field3D &var, Mesh::deriv_func func, CELL_L result = this->fromFieldAligned(result); } - result.setLocation(outloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false; @@ -1156,6 +1152,7 @@ const Field3D Mesh::applyZdiff(const Field3D &var, Mesh::deriv_func func, CELL_L Field3D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); // Check that the input variable has data ASSERT1(var.isAllocated()); @@ -1174,8 +1171,6 @@ const Field3D Mesh::applyZdiff(const Field3D &var, Mesh::deriv_func func, CELL_L } } - result.setLocation(outloc); - return result; } @@ -1197,8 +1192,6 @@ const Field3D Mesh::indexDDX(const Field3D &f, CELL_LOC outloc, DIFF_METHOD meth ASSERT1(outloc == inloc || (outloc == CELL_CENTRE && inloc == CELL_XLOW) || (outloc == CELL_XLOW && inloc == CELL_CENTRE)); - Field3D result(this); - if (this->StaggerGrids && (outloc != inloc)) { // Shifting in X. Centre -> Xlow, or Xlow -> Centre @@ -1213,9 +1206,7 @@ const Field3D Mesh::indexDDX(const Field3D &f, CELL_LOC outloc, DIFF_METHOD meth throw BoutException("Cannot use FFT for X derivatives"); } - result = applyXdiff(f, func, outloc, region); - - return result; + return applyXdiff(f, func, outloc, region); } const Field2D Mesh::indexDDX(const Field2D &f, CELL_LOC outloc, @@ -1239,8 +1230,6 @@ const Field3D Mesh::indexDDY(const Field3D &f, CELL_LOC outloc, DIFF_METHOD meth ASSERT1(outloc == inloc || (outloc == CELL_CENTRE && inloc == CELL_YLOW) || (outloc == CELL_YLOW && inloc == CELL_CENTRE)); - Field3D result(this); - if (this->StaggerGrids && (outloc != inloc)) { // Shifting in Y. Centre -> Ylow, or Ylow -> Centre func = sfDDY; // Set default @@ -1254,9 +1243,7 @@ const Field3D Mesh::indexDDY(const Field3D &f, CELL_LOC outloc, DIFF_METHOD meth throw BoutException("Cannot use FFT for Y derivatives"); } - result = applyYdiff(f, func, outloc, region); - - return result; + return applyYdiff(f, func, outloc, region); } const Field2D Mesh::indexDDY(const Field2D &f, CELL_LOC outloc, @@ -1420,8 +1407,6 @@ const Field3D Mesh::indexD2DX2(const Field3D &f, CELL_LOC outloc, ASSERT1(this == f.getMesh()); - Field3D result(this); - if (StaggerGrids && (outloc != inloc)) { // Shifting in X. Centre -> Xlow, or Xlow -> Centre func = sfD2DX2; // Set default @@ -1435,9 +1420,7 @@ const Field3D Mesh::indexD2DX2(const Field3D &f, CELL_LOC outloc, throw BoutException("Cannot use FFT for X derivatives"); } - result = applyXdiff(f, func, outloc, region); - - return result; + return applyXdiff(f, func, outloc, region); } /*! @@ -1483,8 +1466,6 @@ const Field3D Mesh::indexD2DY2(const Field3D &f, CELL_LOC outloc, ASSERT1(outloc == inloc || (outloc == CELL_CENTRE && inloc == CELL_YLOW) || (outloc == CELL_YLOW && inloc == CELL_CENTRE)); - Field3D result(this); - if (StaggerGrids && (outloc != inloc)) { // Shifting in Y. Centre -> Ylow, or Ylow -> Centre func = sfD2DY2; // Set default @@ -1498,9 +1479,7 @@ const Field3D Mesh::indexD2DY2(const Field3D &f, CELL_LOC outloc, throw BoutException("Cannot use FFT for Y derivatives"); } - result = applyYdiff(f, func, outloc, region); - - return result; + return applyYdiff(f, func, outloc, region); } /*! @@ -1726,6 +1705,7 @@ const Field2D Mesh::indexVDDX(const Field2D &v, const Field2D &f, CELL_LOC outlo Field2D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); if (this->xstart > 1) { // Two or more guard cells @@ -1761,8 +1741,6 @@ const Field2D Mesh::indexVDDX(const Field2D &v, const Field2D &f, CELL_LOC outlo result.bndry_xin = result.bndry_xout = false; #endif - result.setLocation(outloc); - return result; } @@ -1776,6 +1754,7 @@ const Field3D Mesh::indexVDDX(const Field3D &v, const Field3D &f, CELL_LOC outlo ASSERT1(this == v.getMesh()); ASSERT1(this == f.getMesh()); + CELL_LOC vloc = v.getLocation(); CELL_LOC inloc = f.getLocation(); // Input location if (outloc == CELL_DEFAULT) @@ -1787,6 +1766,7 @@ const Field3D Mesh::indexVDDX(const Field3D &v, const Field3D &f, CELL_LOC outlo Field3D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); /// Convert REGION enum to a Region string identifier const auto region_str = REGION_STRING(region); @@ -1926,8 +1906,6 @@ const Field3D Mesh::indexVDDX(const Field3D &v, const Field3D &f, CELL_LOC outlo } } - result.setLocation(outloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false; @@ -1957,10 +1935,10 @@ const Field2D Mesh::indexVDDY(const Field2D &v, const Field2D &f, CELL_LOC outlo Field2D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); if (this->LocalNy == 1){ result=0; - result.setLocation(outloc); return result; } @@ -2102,8 +2080,6 @@ const Field2D Mesh::indexVDDY(const Field2D &v, const Field2D &f, CELL_LOC outlo } } - result.setLocation(outloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false; @@ -2131,10 +2107,10 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo Field3D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); if (this->LocalNy == 1){ result=0; - result.setLocation(outloc); return result; } @@ -2287,8 +2263,6 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo } } - result.setLocation(outloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false; @@ -2318,6 +2292,7 @@ const Field3D Mesh::indexVDDZ(const Field3D &v, const Field3D &f, CELL_LOC outlo Field3D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); /// Convert REGION enum to a Region string identifier const auto region_str = REGION_STRING(region); @@ -2388,8 +2363,6 @@ const Field3D Mesh::indexVDDZ(const Field3D &v, const Field3D &f, CELL_LOC outlo } } - result.setLocation(outloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false; @@ -2427,6 +2400,7 @@ const Field2D Mesh::indexFDDX(const Field2D &v, const Field2D &f, CELL_LOC outlo Field2D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); ASSERT1(this == v.getMesh()); ASSERT1(this == f.getMesh()); @@ -2521,6 +2495,7 @@ const Field3D Mesh::indexFDDX(const Field3D &v, const Field3D &f, CELL_LOC outlo Field3D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); /// Convert REGION enum to a Region string identifier const auto region_str = REGION_STRING(region); @@ -2652,8 +2627,6 @@ const Field3D Mesh::indexFDDX(const Field3D &v, const Field3D &f, CELL_LOC outlo } } - result.setLocation(outloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false; @@ -2691,7 +2664,7 @@ const Field2D Mesh::indexFDDY(const Field2D &v, const Field2D &f, CELL_LOC outlo Field2D result(this); result.allocate(); // Make sure data allocated - result.setLocation(f.getLocation()); + result.setLocation(outloc); /// Convert REGION enum to a Region string identifier const auto region_str = REGION_STRING(region); @@ -2735,8 +2708,6 @@ const Field2D Mesh::indexFDDY(const Field2D &v, const Field2D &f, CELL_LOC outlo } } - result.setLocation(outloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_xin = result.bndry_xout = false; @@ -2788,6 +2759,7 @@ const Field3D Mesh::indexFDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo Field3D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); // If *UseUpDown is true, field "*" has distinct yup and ydown fields which // will be used to calculate a derivative along the magnetic field @@ -2866,8 +2838,6 @@ const Field3D Mesh::indexFDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo result = this->fromFieldAligned(result); } - result.setLocation(outloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false; @@ -2918,6 +2888,7 @@ const Field3D Mesh::indexFDDZ(const Field3D &v, const Field3D &f, CELL_LOC outlo Field3D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); /// Convert REGION enum to a Region string identifier const auto region_str = REGION_STRING(region); @@ -2953,7 +2924,6 @@ const Field3D Mesh::indexFDDZ(const Field3D &v, const Field3D &f, CELL_LOC outlo result[i] = func(vval, fval); } } - result.setLocation(outloc); #if CHECK > 0 // Mark boundaries as invalid