diff --git a/include/bout/index_derivs.hxx b/include/bout/index_derivs.hxx index 28b95f659e..a86ca1c67e 100644 --- a/include/bout/index_derivs.hxx +++ b/include/bout/index_derivs.hxx @@ -123,6 +123,85 @@ public: const metaData meta = func.meta; }; +///////////////////////////////////////////////////////////////////////////////// +/// Following code is for dealing with registering a method/methods for all +/// template combinations, in conjunction with the template_combinations code. +///////////////////////////////////////////////////////////////////////////////// + +struct registerMethod { + template + void operator()(Direction, Stagger, FieldTypeContainer, Method) { + AUTO_TRACE(); + using namespace std::placeholders; + + // Now we want to get the actual field type out of the TypeContainer + // used to pass this around + using FieldType = typename FieldTypeContainer::type; + + Method method{}; + + // Note whilst this should be known at compile time using this directly in the + // template parameters below causes problems for old versions of gcc/libstdc++ + // (tested with 4.8.3) so we currently use a hacky workaround. Once we drop + // support for these versions the branching in the case statement below can be + // removed and we can use nGuard directly in the template statement. + const int nGuards = method.meta.nGuards; + + auto& derivativeRegister = DerivativeStore::getInstance(); + + switch (method.meta.derivType) { + case (DERIV::Standard): + case (DERIV::StandardSecond): + case (DERIV::StandardFourth): { + if (nGuards == 1) { + const auto theFunc = std::bind( + // Method to store in function + &Method::template standard, + // Arguments -- first is hidden this of type-bound, others are placeholders + // for input field, output field, region + method, _1, _2, _3); + derivativeRegister.registerDerivative(theFunc, Direction{}, Stagger{}, method); + } else { + const auto theFunc = std::bind( + // Method to store in function + &Method::template standard, + // Arguments -- first is hidden this of type-bound, others are placeholders + // for input field, output field, region + method, _1, _2, _3); + derivativeRegister.registerDerivative(theFunc, Direction{}, Stagger{}, method); + } + break; + } + case (DERIV::Upwind): + case (DERIV::Flux): { + if (nGuards == 1) { + const auto theFunc = std::bind( + // Method to store in function + &Method::template upwindOrFlux, + // Arguments -- first is hidden this of type-bound, others are placeholders + // for input field, output field, region + method, _1, _2, _3, _4); + derivativeRegister.registerDerivative(theFunc, Direction{}, Stagger{}, method); + } else { + const auto theFunc = std::bind( + // Method to store in function + &Method::template upwindOrFlux, + // Arguments -- first is hidden this of type-bound, others are placeholders + // for input field, output field, region + method, _1, _2, _3, _4); + derivativeRegister.registerDerivative(theFunc, Direction{}, Stagger{}, method); + } + break; + } + default: + throw BoutException("Unhandled derivative method in registerMethod."); + }; + } +}; + #define DEFINE_STANDARD_DERIV_CORE(name, key, nGuards, type) \ struct name { \ BoutReal operator()(const stencil& f) const; \ @@ -171,18 +250,82 @@ public: #define DEFINE_FLUX_DERIV_STAGGERED(name, key, nGuards, type) \ DEFINE_FLUX_DERIV(name, key, nGuards, type) +/// Some helper defines for now that allow us to wrap up enums +/// and the specific methods. +#define WRAP_ENUM(family, value) enumWrapper + +#define REGISTER_DERIVATIVE(name) \ + namespace { \ + produceCombinations, \ + Set, \ + Set, TypeContainer>, \ + Set>> \ + reg##name(registerMethod{}); \ + } +#define REGISTER_STAGGERED_DERIVATIVE(name) \ + namespace { \ + produceCombinations, \ + Set, \ + Set, TypeContainer>, \ + Set>> \ + reg##name(registerMethod{}); \ + } + +#define REGISTER_STANDARD_DERIVATIVE(name, key, nGuards, type) \ + DEFINE_STANDARD_DERIV_CORE(name, key, nGuards, type) \ + REGISTER_DERIVATIVE(name) \ + BoutReal name::operator()(const stencil& f) const + +#define REGISTER_UPWIND_DERIVATIVE(name, key, nGuards, type) \ + DEFINE_UPWIND_DERIV_CORE(name, key, nGuards, type) \ + REGISTER_DERIVATIVE(name) \ + BoutReal name::operator()(BoutReal vc, const stencil& f) const + +#define REGISTER_FLUX_DERIVATIVE(name, key, nGuards, type) \ + DEFINE_FLUX_DERIV_CORE(name, key, nGuards, type) \ + REGISTER_DERIVATIVE(name) \ + BoutReal name::operator()(const stencil& v, const stencil& f) const + +#define REGISTER_STANDARD_STAGGERED_DERIVATIVE(name, key, nGuards, type) \ + DEFINE_STANDARD_DERIV_CORE(name, key, nGuards, type) \ + REGISTER_STAGGERED_DERIVATIVE(name) \ + BoutReal name::operator()(const stencil& f) const + +#define REGISTER_STANDARD_DERIVATIVE_STAGGERED(name, key, nGuards, type) \ + REGISTER_STANDARD_STAGGERED_DERIVATIVE(name, key, nGuards, type) + +#define REGISTER_UPWIND_STAGGERED_DERIVATIVE(name, key, nGuards, type) \ + /*Note staggered upwind looks like flux*/ \ + DEFINE_FLUX_DERIV_CORE(name, key, nGuards, type) \ + REGISTER_STAGGERED_DERIVATIVE(name) \ + BoutReal name::operator()(const stencil& v, const stencil& f) const + +#define REGISTER_UPWIND_DERIVATIVE_STAGGERED(name, key, nGuards, type) \ + REGISTER_UPWIND_STAGGERED_DERIVATIVE(name, key, nGuards, type) + +#define REGISTER_FLUX_STAGGERED_DERIVATIVE(name, key, nGuards, type) \ + DEFINE_FLUX_DERIV_CORE(name, key, nGuards, type) \ + REGISTER_STAGGERED_DERIVATIVE(name) \ + BoutReal name::operator()(const stencil& v, const stencil& f) const + +#define REGISTER_FLUX_DERIVATIVE_STAGGERED(name, key, nGuards, type) \ + REGISTER_FLUX_STAGGERED_DERIVATIVE(name, key, nGuards, type) ////////////////////// FIRST DERIVATIVES ///////////////////// /// central, 2nd order -DEFINE_STANDARD_DERIV(DDX_C2, "C2", 1, DERIV::Standard) { return 0.5 * (f.p - f.m); }; +REGISTER_STANDARD_DERIVATIVE(DDX_C2, "C2", 1, DERIV::Standard) { + return 0.5 * (f.p - f.m); +}; /// central, 4th order -DEFINE_STANDARD_DERIV(DDX_C4, "C4", 2, DERIV::Standard) { +REGISTER_STANDARD_DERIVATIVE(DDX_C4, "C4", 2, DERIV::Standard) { return (8. * f.p - 8. * f.m + f.mm - f.pp) / 12.; } /// Central WENO method, 2nd order (reverts to 1st order near shocks) -DEFINE_STANDARD_DERIV(DDX_CWENO2, "W2", 1, DERIV::Standard) { +REGISTER_STANDARD_DERIVATIVE(DDX_CWENO2, "W2", 1, DERIV::Standard) { BoutReal isl, isr, isc; // Smoothness indicators BoutReal al, ar, ac, sa; // Un-normalised weights BoutReal dl, dr, dc; // Derivatives using different stencils @@ -204,7 +347,7 @@ DEFINE_STANDARD_DERIV(DDX_CWENO2, "W2", 1, DERIV::Standard) { } // Smoothing 2nd order derivative -DEFINE_STANDARD_DERIV(DDX_S2, "S2", 2, DERIV::Standard) { +REGISTER_STANDARD_DERIVATIVE(DDX_S2, "S2", 2, DERIV::Standard) { // 4th-order differencing BoutReal result = (8. * f.p - 8. * f.m + f.mm - f.pp) / 12.; @@ -221,19 +364,19 @@ DEFINE_STANDARD_DERIV(DDX_S2, "S2", 2, DERIV::Standard) { ////////////////////////////// /// Second derivative: Central, 2nd order -DEFINE_STANDARD_DERIV(D2DX2_C2, "C2", 1, DERIV::StandardSecond) { +REGISTER_STANDARD_DERIVATIVE(D2DX2_C2, "C2", 1, DERIV::StandardSecond) { return f.p + f.m - 2. * f.c; } /// Second derivative: Central, 4th order -DEFINE_STANDARD_DERIV(D2DX2_C4, "C4", 2, DERIV::StandardSecond) { +REGISTER_STANDARD_DERIVATIVE(D2DX2_C4, "C4", 2, DERIV::StandardSecond) { return (-f.pp + 16. * f.p - 30. * f.c + 16. * f.m - f.mm) / 12.; } ////////////////////////////// //--- Fourth order derivatives ////////////////////////////// -DEFINE_STANDARD_DERIV(D4DX4_C2, "C2", 2, DERIV::StandardFourth) { +REGISTER_STANDARD_DERIVATIVE(D4DX4_C2, "C2", 2, DERIV::StandardFourth) { return (f.pp - 4. * f.p + 6. * f.c - 4. * f.m + f.mm); } @@ -250,15 +393,17 @@ std::tuple vUpDown(BoutReal v) { } /// Upwinding: Central, 2nd order -DEFINE_UPWIND_DERIV(VDDX_C2, "C2", 1, DERIV::Upwind) { return vc * 0.5 * (f.p - f.m); } +REGISTER_UPWIND_DERIVATIVE(VDDX_C2, "C2", 1, DERIV::Upwind) { + return vc * 0.5 * (f.p - f.m); +} /// Upwinding: Central, 4th order -DEFINE_UPWIND_DERIV(VDDX_C4, "C4", 2, DERIV::Upwind) { +REGISTER_UPWIND_DERIVATIVE(VDDX_C4, "C4", 2, DERIV::Upwind) { return vc * (8. * f.p - 8. * f.m + f.mm - f.pp) / 12.; } /// upwind, 1st order -DEFINE_UPWIND_DERIV(VDDX_U1, "U1", 1, DERIV::Upwind) { // No vec +REGISTER_UPWIND_DERIVATIVE(VDDX_U1, "U1", 1, DERIV::Upwind) { // No vec // Existing form doesn't vectorise due to branching return vc >= 0.0 ? vc * (f.c - f.m) : vc * (f.p - f.c); // Alternative form would but may involve more operations @@ -267,7 +412,7 @@ DEFINE_UPWIND_DERIV(VDDX_U1, "U1", 1, DERIV::Upwind) { // No vec } /// upwind, 2nd order -DEFINE_UPWIND_DERIV(VDDX_U2, "U2", 2, DERIV::Upwind) { // No vec +REGISTER_UPWIND_DERIVATIVE(VDDX_U2, "U2", 2, DERIV::Upwind) { // No vec // Existing form doesn't vectorise due to branching return vc >= 0.0 ? vc * (1.5 * f.c - 2.0 * f.m + 0.5 * f.mm) : vc * (-0.5 * f.pp + 2.0 * f.p - 1.5 * f.c); @@ -278,7 +423,7 @@ DEFINE_UPWIND_DERIV(VDDX_U2, "U2", 2, DERIV::Upwind) { // No vec } /// upwind, 3rd order -DEFINE_UPWIND_DERIV(VDDX_U3, "U3", 2, DERIV::Upwind) { // No vec +REGISTER_UPWIND_DERIVATIVE(VDDX_U3, "U3", 2, DERIV::Upwind) { // No vec // Existing form doesn't vectorise due to branching return vc >= 0.0 ? vc * (4. * f.p - 12. * f.m + 2. * f.mm + 6. * f.c) / 12. : vc * (-4. * f.m + 12. * f.p - 2. * f.pp - 6. * f.c) / 12.; @@ -290,7 +435,7 @@ DEFINE_UPWIND_DERIV(VDDX_U3, "U3", 2, DERIV::Upwind) { // No vec } /// 3rd-order WENO scheme -DEFINE_UPWIND_DERIV(VDDX_WENO3, "W3", 2, DERIV::Upwind) { // No vec +REGISTER_UPWIND_DERIVATIVE(VDDX_WENO3, "W3", 2, DERIV::Upwind) { // No vec BoutReal deriv, w, r; // Existing form doesn't vectorise due to branching @@ -319,7 +464,7 @@ DEFINE_UPWIND_DERIV(VDDX_WENO3, "W3", 2, DERIV::Upwind) { // No vec ///----------------------------------------------------------------- /// 3rd-order CWENO. Uses the upwinding code and split flux -DEFINE_STANDARD_DERIV(DDX_CWENO3, "W3", 2, DERIV::Standard) { +REGISTER_STANDARD_DERIVATIVE(DDX_CWENO3, "W3", 2, DERIV::Standard) { BoutReal a, ma = fabs(f.c); // Split flux a = fabs(f.m); @@ -362,7 +507,7 @@ DEFINE_STANDARD_DERIV(DDX_CWENO3, "W3", 2, DERIV::Standard) { /// //////////////////////////////////////////////////////////////////////////////// -DEFINE_FLUX_DERIV(FDDX_U1, "U1", 1, DERIV::Flux) { // No vec +REGISTER_FLUX_DERIVATIVE(FDDX_U1, "U1", 1, DERIV::Flux) { // No vec // Velocity at lower end BoutReal vs = 0.5 * (v.m + v.c); @@ -378,9 +523,11 @@ DEFINE_FLUX_DERIV(FDDX_U1, "U1", 1, DERIV::Flux) { // No vec return result - std::get<0>(vSplit) * f.c + std::get<1>(vSplit) * f.p; } -DEFINE_FLUX_DERIV(FDDX_C2, "C2", 2, DERIV::Flux) { return 0.5 * (v.p * f.p - v.m * f.m); } +REGISTER_FLUX_DERIVATIVE(FDDX_C2, "C2", 2, DERIV::Flux) { + return 0.5 * (v.p * f.p - v.m * f.m); +} -DEFINE_FLUX_DERIV(FDDX_C4, "C4", 2, DERIV::Flux) { +REGISTER_FLUX_DERIVATIVE(FDDX_C4, "C4", 2, DERIV::Flux) { return (8. * v.p * f.p - 8. * v.m * f.m + v.mm * f.mm - v.pp * f.pp) / 12.; } @@ -407,25 +554,25 @@ DEFINE_FLUX_DERIV(FDDX_C4, "C4", 2, DERIV::Flux) { //////////////////////////////////////////////////////////////////////////////// /// Standard methods -- first order //////////////////////////////////////////////////////////////////////////////// -DEFINE_STANDARD_DERIV_STAGGERED(DDX_C2_stag, "C2", 1, DERIV::Standard) { +REGISTER_STANDARD_DERIVATIVE_STAGGERED(DDX_C2_stag, "C2", 1, DERIV::Standard) { return f.p - f.m; } -DEFINE_STANDARD_DERIV_STAGGERED(DDX_C4_stag, "C4", 2, DERIV::Standard) { +REGISTER_STANDARD_DERIVATIVE_STAGGERED(DDX_C4_stag, "C4", 2, DERIV::Standard) { return (27. * (f.p - f.m) - (f.pp - f.mm)) / 24.; } //////////////////////////////////////////////////////////////////////////////// /// Standard methods -- second order //////////////////////////////////////////////////////////////////////////////// -DEFINE_STANDARD_DERIV_STAGGERED(D2DX2_C2_stag, "C2", 2, DERIV::StandardSecond) { +REGISTER_STANDARD_DERIVATIVE_STAGGERED(D2DX2_C2_stag, "C2", 2, DERIV::StandardSecond) { return (f.pp + f.mm - f.p - f.m) / 2.; } //////////////////////////////////////////////////////////////////////////////// /// Upwind methods //////////////////////////////////////////////////////////////////////////////// -DEFINE_UPWIND_DERIV_STAGGERED(VDDX_U1_stag, "U1", 1, DERIV::Upwind) { +REGISTER_UPWIND_DERIVATIVE_STAGGERED(VDDX_U1_stag, "U1", 1, DERIV::Upwind) { // Lower cell boundary BoutReal result = (v.m >= 0) ? v.m * f.m : v.m * f.c; @@ -438,7 +585,7 @@ DEFINE_UPWIND_DERIV_STAGGERED(VDDX_U1_stag, "U1", 1, DERIV::Upwind) { return result; } -DEFINE_UPWIND_DERIV_STAGGERED(VDDX_U2_stag, "U2", 2, DERIV::Upwind) { +REGISTER_UPWIND_DERIVATIVE_STAGGERED(VDDX_U2_stag, "U2", 2, DERIV::Upwind) { // Calculate d(v*f)/dx = (v*f)[i+1/2] - (v*f)[i-1/2] // Upper cell boundary @@ -454,13 +601,13 @@ DEFINE_UPWIND_DERIV_STAGGERED(VDDX_U2_stag, "U2", 2, DERIV::Upwind) { return result; } -DEFINE_UPWIND_DERIV_STAGGERED(VDDX_C2_stag, "C2", 1, DERIV::Upwind) { +REGISTER_UPWIND_DERIVATIVE_STAGGERED(VDDX_C2_stag, "C2", 1, DERIV::Upwind) { // Result is needed at location of f: interpolate v to f's location and take an // unstaggered derivative of f return 0.5 * (v.p + v.m) * 0.5 * (f.p - f.m); } -DEFINE_UPWIND_DERIV_STAGGERED(VDDX_C4_stag, "C4", 2, DERIV::Upwind) { +REGISTER_UPWIND_DERIVATIVE_STAGGERED(VDDX_C4_stag, "C4", 2, DERIV::Upwind) { // Result is needed at location of f: interpolate v to f's location and take an // unstaggered derivative of f return (9. * (v.m + v.p) - v.mm - v.pp) / 16. * (8. * f.p - 8. * f.m + f.mm - f.pp) @@ -470,7 +617,7 @@ DEFINE_UPWIND_DERIV_STAGGERED(VDDX_C4_stag, "C4", 2, DERIV::Upwind) { //////////////////////////////////////////////////////////////////////////////// /// Flux methods //////////////////////////////////////////////////////////////////////////////// -DEFINE_FLUX_DERIV_STAGGERED(FDDX_U1_stag, "U1", 1, DERIV::Flux) { +REGISTER_FLUX_DERIVATIVE_STAGGERED(FDDX_U1_stag, "U1", 1, DERIV::Flux) { // Lower cell boundary BoutReal result = (v.m >= 0) ? v.m * f.m : v.m * f.c; @@ -480,208 +627,10 @@ DEFINE_FLUX_DERIV_STAGGERED(FDDX_U1_stag, "U1", 1, DERIV::Flux) { return -result; } -///////////////////////////////////////////////////////////////////////////////// -/// Following code is for dealing with registering a method/methods for all -/// template combinations, in conjunction with the template_combinations code. -///////////////////////////////////////////////////////////////////////////////// - -struct registerMethod { - template - void operator()(Direction, Stagger, FieldTypeContainer, Method) { - AUTO_TRACE(); - using namespace std::placeholders; - - // Now we want to get the actual field type out of the TypeContainer - // used to pass this around - using FieldType = typename FieldTypeContainer::type; - - Method method{}; - - // Note whilst this should be known at compile time using this directly in the - // template parameters below causes problems for old versions of gcc/libstdc++ - // (tested with 4.8.3) so we currently use a hacky workaround. Once we drop - // support for these versions the branching in the case statement below can be - // removed and we can use nGuard directly in the template statement. - const int nGuards = method.meta.nGuards; - - auto& derivativeRegister = DerivativeStore::getInstance(); - - switch (method.meta.derivType) { - case (DERIV::Standard): - case (DERIV::StandardSecond): - case (DERIV::StandardFourth): { - if (nGuards == 1) { - const auto theFunc = std::bind( - // Method to store in function - &Method::template standard, - // Arguments -- first is hidden this of type-bound, others are placeholders - // for input field, output field, region - method, _1, _2, _3); - derivativeRegister.registerDerivative(theFunc, Direction{}, Stagger{}, method); - } else { - const auto theFunc = std::bind( - // Method to store in function - &Method::template standard, - // Arguments -- first is hidden this of type-bound, others are placeholders - // for input field, output field, region - method, _1, _2, _3); - derivativeRegister.registerDerivative(theFunc, Direction{}, Stagger{}, method); - } - break; - } - case (DERIV::Upwind): - case (DERIV::Flux): { - if (nGuards == 1) { - const auto theFunc = std::bind( - // Method to store in function - &Method::template upwindOrFlux, - // Arguments -- first is hidden this of type-bound, others are placeholders - // for input field, output field, region - method, _1, _2, _3, _4); - derivativeRegister.registerDerivative(theFunc, Direction{}, Stagger{}, method); - } else { - const auto theFunc = std::bind( - // Method to store in function - &Method::template upwindOrFlux, - // Arguments -- first is hidden this of type-bound, others are placeholders - // for input field, output field, region - method, _1, _2, _3, _4); - derivativeRegister.registerDerivative(theFunc, Direction{}, Stagger{}, method); - } - break; - } - default: - throw BoutException("Unhandled derivative method in registerMethod."); - }; - } -}; - -/// Some helper defines for now that allow us to wrap up enums -/// and the specific methods. -#define WRAP_ENUM(family, value) enumWrapper - -#define REGISTER_DERIVATIVE(name) \ - namespace { \ - produceCombinations, \ - Set, \ - Set, TypeContainer>, \ - Set>> \ - reg(registerMethod{}); \ - } -#define REGISTER_STAGGERED_DERIVATIVE(name) \ - namespace { \ - produceCombinations, \ - Set, \ - Set, TypeContainer>, \ - Set>> \ - reg(registerMethod{}); \ - } - -#define REGISTER_STANDARD_DERIVATIVE(name, key, nGuards, type) \ - DEFINE_STANDARD_DERIV_CORE(name, key, nGuards, type) \ - REGISTER_DERIVATIVE(name) \ - BoutReal name::operator()(const stencil& f) const - -#define REGISTER_UPWIND_DERIVATIVE(name, key, nGuards, type) \ - DEFINE_UPWIND_DERIV_CORE(name, key, nGuards, type) \ - REGISTER_DERIVATIVE(name) \ - BoutReal name::operator()(BoutReal vc, const stencil& f) const - -#define REGISTER_FLUX_DERIVATIVE(name, key, nGuards, type) \ - DEFINE_FLUX_DERIV_CORE(name, key, nGuards, type) \ - REGISTER_DERIVATIVE(name) \ - BoutReal name::operator()(const stencil& v, const stencil& f) const - -#define REGISTER_STANDARD_STAGGERED_DERIVATIVE(name, key, nGuards, type) \ - DEFINE_STANDARD_DERIV_CORE(name, key, nGuards, type) \ - REGISTER_STAGGERED_DERIVATIVE(name) \ - BoutReal name::operator()(const stencil& f) const - -#define REGISTER_UPWIND_STAGGERED_DERIVATIVE(name, key, nGuards, type) \ - /*Note staggered upwind looks like flux*/ \ - DEFINE_FLUX_DERIV_CORE(name, key, nGuards, type) \ - REGISTER_STAGGERED_DERIVATIVE(name) \ - BoutReal name::operator()(const stencil& v, const stencil& f) const - -#define REGISTER_FLUX_STAGGERED_DERIVATIVE(name, key, nGuards, type) \ - DEFINE_FLUX_DERIV_CORE(name, key, nGuards, type) \ - REGISTER_STAGGERED_DERIVATIVE(name) \ - BoutReal name::operator()(const stencil& v, const stencil& f) const - -///////////////////////////////////////////////////////////////////////////////// -/// Here's an example of registering a couple of DerivativeType methods -/// at once for no staggering -///////////////////////////////////////////////////////////////////////////////// - -produceCombinations, - Set, - Set, TypeContainer>, - Set< - // Standard - DerivativeType, DerivativeType, - DerivativeType, DerivativeType, - DerivativeType, - // Standard 2nd order - DerivativeType, DerivativeType, - // Standard 4th order - DerivativeType, - // Upwind - DerivativeType, DerivativeType, - DerivativeType, DerivativeType, - DerivativeType, DerivativeType, - // Flux - DerivativeType, DerivativeType, - DerivativeType>> - registerDerivatives(registerMethod{}); - -produceCombinations, Set, - Set, TypeContainer>, - Set< - // Standard - DerivativeType, DerivativeType, - // Standard 2nd order - DerivativeType, - // Standard 4th order - // Upwind - DerivativeType, DerivativeType, - // Flux - DerivativeType>> - registerDerivativesYOrtho(registerMethod{}); - -produceCombinations, - Set, - Set, TypeContainer>, - Set< - // Standard - DerivativeType, DerivativeType, - // Standard 2nd order - DerivativeType, - // Upwind - DerivativeType, DerivativeType, - DerivativeType, DerivativeType, - // Flux - DerivativeType>> - registerStaggeredDerivatives(registerMethod{}); - -produceCombinations, - Set, - Set, TypeContainer>, - Set< - // Standard - DerivativeType, - // Standard 2nd order - // Upwind - DerivativeType, DerivativeType, - // Flux - DerivativeType>> - registerStaggeredDerivativesYOrtho(registerMethod{}); +///////////////////////////////////////////////////////////////////////////////////// +/// Here's an example of defining and registering a custom method that doesn't fit +/// into the standard stencil based approach. +// ///////////////////////////////////////////////////////////////////////////////// class FFTDerivativeType { public: diff --git a/include/bout/index_derivs_interface.hxx b/include/bout/index_derivs_interface.hxx index 154a6b9fd6..4ebdc2ac59 100644 --- a/include/bout/index_derivs_interface.hxx +++ b/include/bout/index_derivs_interface.hxx @@ -207,14 +207,13 @@ template T DDY(const T& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT", REGION region = RGN_NOBNDRY) { AUTO_TRACE(); - if (std::is_base_of::value && f.hasYupYdown() - && ((&f.yup() != &f) || (&f.ydown() != &f))) { + if (std::is_base_of::value && f.hasYupYdown()) { return standardDerivative(f, outloc, method, region); } else { const T f_aligned = f.getMesh()->toFieldAligned(f); - T result = - standardDerivative(f_aligned, outloc, method, region); + T result = standardDerivative(f_aligned, outloc, + method, region); return f.getMesh()->fromFieldAligned(result); } } @@ -223,14 +222,13 @@ template T D2DY2(const T& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT", REGION region = RGN_NOBNDRY) { AUTO_TRACE(); - if (std::is_base_of::value && f.hasYupYdown() - && ((&f.yup() != &f) || (&f.ydown() != &f))) { + if (std::is_base_of::value && f.hasYupYdown()) { return standardDerivative( f, outloc, method, region); } else { const T f_aligned = f.getMesh()->toFieldAligned(f); - T result = standardDerivative(f_aligned, outloc, - method, region); + T result = standardDerivative( + f_aligned, outloc, method, region); return f.getMesh()->fromFieldAligned(result); } } @@ -239,14 +237,13 @@ template T D4DY4(const T& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT", REGION region = RGN_NOBNDRY) { AUTO_TRACE(); - if (std::is_base_of::value && f.hasYupYdown() - && ((&f.yup() != &f) || (&f.ydown() != &f))) { + if (std::is_base_of::value && f.hasYupYdown()) { return standardDerivative( f, outloc, method, region); } else { const T f_aligned = f.getMesh()->toFieldAligned(f); - T result = standardDerivative(f_aligned, outloc, - method, region); + T result = standardDerivative( + f_aligned, outloc, method, region); return f.getMesh()->fromFieldAligned(result); } } @@ -312,10 +309,8 @@ template T VDDY(const T& vel, const T& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT", REGION region = RGN_NOBNDRY) { AUTO_TRACE(); - bool fHasParallelSlices = (std::is_base_of::value && f.hasYupYdown() - && ((&f.yup() != &f) || (&f.ydown() != &f))); - bool velHasParallelSlices = (std::is_base_of::value && vel.hasYupYdown() - && ((&vel.yup() != &vel) || (&vel.ydown() != &vel))); + bool fHasParallelSlices = (std::is_base_of::value && f.hasYupYdown()); + bool velHasParallelSlices = (std::is_base_of::value && vel.hasYupYdown()); if (fHasParallelSlices && velHasParallelSlices) { return flowDerivative(vel, f, outloc, method, region); @@ -332,10 +327,8 @@ template T FDDY(const T& vel, const T& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT", REGION region = RGN_NOBNDRY) { AUTO_TRACE(); - bool fHasParallelSlices = (std::is_base_of::value && f.hasYupYdown() - && ((&f.yup() != &f) || (&f.ydown() != &f))); - bool velHasParallelSlices = (std::is_base_of::value && vel.hasYupYdown() - && ((&vel.yup() != &vel) || (&vel.ydown() != &vel))); + bool fHasParallelSlices = (std::is_base_of::value && f.hasYupYdown()); + bool velHasParallelSlices = (std::is_base_of::value && vel.hasYupYdown()); if (fHasParallelSlices && velHasParallelSlices) { return flowDerivative(vel, f, outloc, method, region); diff --git a/include/bout/paralleltransform.hxx b/include/bout/paralleltransform.hxx index 2f85136be4..c7b6ebf7d8 100644 --- a/include/bout/paralleltransform.hxx +++ b/include/bout/paralleltransform.hxx @@ -59,7 +59,7 @@ public: * Merges the yup and ydown() fields of f, so that * f.yup() = f.ydown() = f */ - void calcYUpDown(Field3D &f) override {f.mergeYupYdown();} + void calcYUpDown(Field3D &f) override; /*! * The field is already aligned in Y, so this @@ -125,7 +125,7 @@ public: } /// A 3D array, implemented as nested vectors - using arr3Dvec = std::vector>>; + using arr3Dvec = std::vector>>; private: Mesh &mesh; ///< The mesh this paralleltransform is part of @@ -137,8 +137,19 @@ private: /// Cache of phase shifts for transforming from field-aligned coordinates to X-Z orthogonal coordinates arr3Dvec fromAlignedPhs; - arr3Dvec yupPhs; ///< Cache of phase shifts for calculating yup fields - arr3Dvec ydownPhs; ///< Cache of phase shifts for calculating ydown fields + /// Helper POD for parallel slice phase shifts + struct ParallelSlicePhase { + arr3Dvec phase_shift; + int y_offset; + }; + + /// Cache of phase shifts for the parallel slices. Slices are stored + /// in the following order: + /// {+1, ..., +n, -1, ..., -n} + /// slice[i] stores offset i+1 + /// slice[2*i + 1] stores offset -(i+1) + /// where i goes from 0 to (n-1), with n the number of y guard cells + std::vector parallel_slice_phases; /*! * Shift a 2D field in Z. @@ -185,11 +196,19 @@ private: * @param[in] phs Phase shift, assumed to have length (mesh.LocalNz/2 + 1) i.e. the number of modes * @param[out] out A 1D array of length mesh.LocalNz, already allocated */ - void shiftZ(const BoutReal *in, const std::vector &phs, BoutReal *out) const; + void shiftZ(const BoutReal *in, const Array &phs, BoutReal *out) const; /// Calculate and store the phases for to/from field aligned and for /// the parallel slices using zShift void cachePhases(); + + /// Shift a 3D field \p f in Z to all the parallel slices in \p phases + /// + /// @param[in] f The field to shift + /// @param[in] phases The phase and offset information for each parallel slice + /// @return The shifted parallel slices + std::vector shiftZ(const Field3D& f, + const std::vector& phases) const; }; diff --git a/include/field2d.hxx b/include/field2d.hxx index adbf8b0800..eae5befe62 100644 --- a/include/field2d.hxx +++ b/include/field2d.hxx @@ -129,7 +129,10 @@ class Field2D : public Field, public FieldData { const Field2D& ydown() const { return *this; } - + + Field2D& ynext(int UNUSED(dir)) { return *this; } + const Field2D& ynext(int UNUSED(dir)) const { return *this; } + // Operators /*! diff --git a/include/field3d.hxx b/include/field3d.hxx index f72d21f502..f9b23400ba 100644 --- a/include/field3d.hxx +++ b/include/field3d.hxx @@ -40,6 +40,8 @@ class Mesh; // #include "bout/mesh.hxx" #include "bout/field_visitor.hxx" +#include + /// Class for 3D X-Y-Z scalar fields /*! This class represents a scalar field defined over the mesh. @@ -228,35 +230,37 @@ class Field3D : public Field, public FieldData { /// Check if this field has yup and ydown fields bool hasYupYdown() const { - return (yup_field != nullptr) && (ydown_field != nullptr); + return !yup_fields.empty() and !ydown_fields.empty(); } /// Return reference to yup field - Field3D& yup() { - ASSERT2(yup_field != nullptr); // Check for communicate - return *yup_field; + Field3D &yup(std::vector::size_type index = 0) { + ASSERT2(index < yup_fields.size()); + return yup_fields[index]; } /// Return const reference to yup field - const Field3D& yup() const { - ASSERT2(yup_field != nullptr); - return *yup_field; + const Field3D &yup(std::vector::size_type index = 0) const { + ASSERT2(index < yup_fields.size()); + return yup_fields[index]; } - + /// Return reference to ydown field - Field3D& ydown() { - ASSERT2(ydown_field != nullptr); - return *ydown_field; + Field3D &ydown(std::vector::size_type index = 0) { + ASSERT2(index < ydown_fields.size()); + return ydown_fields[index]; } - + /// Return const reference to ydown field - const Field3D& ydown() const { - ASSERT2(ydown_field != nullptr); - return *ydown_field; + const Field3D &ydown(std::vector::size_type index = 0) const { + ASSERT2(index < ydown_fields.size()); + return ydown_fields[index]; } - /// Return yup if dir=+1, and ydown if dir=-1 - Field3D& ynext(int dir); - const Field3D& ynext(int dir) const; + /// Return the parallel slice at \p offset + /// + /// \p offset of 0 returns the main field itself + Field3D& ynext(int offset); + const Field3D& ynext(int offset) const; /// Set variable location for staggered grids to @param new_location /// @@ -462,8 +466,8 @@ class Field3D : public Field, public FieldData { swap(first.nz, second.nz); swap(first.location, second.location); swap(first.deriv, second.deriv); - swap(first.yup_field, second.yup_field); - swap(first.ydown_field, second.ydown_field); + swap(first.yup_fields, second.yup_fields); + swap(first.ydown_fields, second.ydown_fields); swap(first.bndry_op, second.bndry_op); swap(first.boundaryIsCopy, second.boundaryIsCopy); swap(first.boundaryIsSet, second.boundaryIsSet); @@ -487,8 +491,8 @@ private: /// Time derivative (may be nullptr) Field3D *deriv{nullptr}; - /// Pointers to fields containing values along Y - Field3D *yup_field{nullptr}, *ydown_field{nullptr}; + /// Fields containing values along Y + std::vector yup_fields{}, ydown_fields{}; }; // Non-member overloaded operators diff --git a/include/interpolation.hxx b/include/interpolation.hxx index 6c688afb69..b20389a00e 100644 --- a/include/interpolation.hxx +++ b/include/interpolation.hxx @@ -118,61 +118,38 @@ const T interp_to(const T& var, CELL_LOC loc, REGION region = RGN_ALL) { // At least 2 boundary cells needed for interpolation in y-direction ASSERT0(fieldmesh->ystart >= 2); - if (var.hasYupYdown() && ((&var.yup() != &var) || (&var.ydown() != &var))) { - // Field "var" has distinct yup and ydown fields which - // will be used to calculate a derivative along - // the magnetic field - throw BoutException( - "At the moment, fields with yup/ydown cannot use interp_to.\n" - "If we implement a 3-point stencil for interpolate or double-up\n" - "/double-down fields, then we can use this case."); - - if ((location == CELL_CENTRE) && (loc == CELL_YLOW)) { // C2L - BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { - // Producing a stencil centred around a lower X value - result[i] = interp( - populateStencil(var, i)); - } - } else if (location == CELL_YLOW) { // L2C - BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { - // Stencil centred around a cell centre - result[i] = interp( - populateStencil(var, i)); - } + // We can't interpolate in y unless we're field-aligned + // FIXME: Add check once we label fields as orthogonal/aligned + + const T var_fa = fieldmesh->toFieldAligned(var); + if (region != RGN_NOBNDRY) { + // repeat the hack above for boundary points + // this avoids a duplicate toFieldAligned call if we had called + // result = toFieldAligned(result) + // to get the boundary cells + // + // result is requested in some boundary region(s) + result = var_fa; // NOTE: This is just for boundaries. FIX! + result.allocate(); + result.setLocation(loc); + } + + if ((location == CELL_CENTRE) && (loc == CELL_YLOW)) { // C2L + BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { + // Producing a stencil centred around a lower X value + result[i] = + interp(populateStencil(var_fa, i)); } - } else { - // var has no yup/ydown fields, so we need to shift into field-aligned - // coordinates - - const T var_fa = fieldmesh->toFieldAligned(var); - if (region != RGN_NOBNDRY) { - // repeat the hack above for boundary points - // this avoids a duplicate toFieldAligned call if we had called - // result = toFieldAligned(result) - // to get the boundary cells - // - // result is requested in some boundary region(s) - result = var_fa; // NOTE: This is just for boundaries. FIX! - result.allocate(); - result.setLocation(loc); + } else if (location == CELL_YLOW) { // L2C + BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { + // Stencil centred around a cell centre + result[i] = + interp(populateStencil(var_fa, i)); } + } - if ((location == CELL_CENTRE) && (loc == CELL_YLOW)) { // C2L - BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { - // Producing a stencil centred around a lower X value - result[i] = interp( - populateStencil(var_fa, i)); - } - } else if (location == CELL_YLOW) { // L2C - BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { - // Stencil centred around a cell centre - result[i] = interp( - populateStencil(var_fa, i)); - } - } + result = fieldmesh->fromFieldAligned(result); - result = fieldmesh->fromFieldAligned(result); - } break; } case CELL_ZLOW: { diff --git a/include/stencils.hxx b/include/stencils.hxx index 517929b1dc..fe3b8b1160 100644 --- a/include/stencils.hxx +++ b/include/stencils.hxx @@ -46,49 +46,73 @@ void inline populateStencil(stencil &s, const FieldType& f, const typename Field switch(stagger) { case(STAGGER::None): - if (nGuard == 2) s.mm = f[i.template minus<2, direction>()]; + if (nGuard == 2) { + if (direction == DIRECTION::YOrthogonal) { + s.mm = f.ynext(-2)[i.template minus<2, direction>()]; + } else { + s.mm = f[i.template minus<2, direction>()]; + } + } if (direction == DIRECTION::YOrthogonal) { - s.m = f.ydown()[i.template minus<1, direction>()]; + s.m = f.ynext(-1)[i.template minus<1, direction>()]; } else { s.m = f[i.template minus<1, direction>()]; } s.c = f[i]; if (direction == DIRECTION::YOrthogonal) { - s.p = f.yup()[i.template plus<1, direction>()]; + s.p = f.ynext(1)[i.template plus<1, direction>()]; } else { s.p = f[i.template plus<1, direction>()]; } - if (nGuard == 2) s.pp = f[i.template plus<2, direction>()]; + if (nGuard == 2) { + if (direction == DIRECTION::YOrthogonal) { + s.pp = f.ynext(2)[i.template plus<2, direction>()]; + } else { + s.pp = f[i.template plus<2, direction>()]; + } + } break; case(STAGGER::C2L): - if (nGuard == 2) s.mm = f[i.template minus<2, direction>()]; + if (nGuard == 2) { + if (direction == DIRECTION::YOrthogonal) { + s.mm = f.ynext(-2)[i.template minus<2, direction>()]; + } else { + s.mm = f[i.template minus<2, direction>()]; + } + } if (direction == DIRECTION::YOrthogonal) { - s.m = f.ydown()[i.template minus<1, direction>()]; + s.m = f.ynext(-1)[i.template minus<1, direction>()]; } else { s.m = f[i.template minus<1, direction>()]; } s.c = f[i]; s.p = s.c; if (direction == DIRECTION::YOrthogonal) { - s.pp = f.yup()[i.template plus<1, direction>()]; + s.pp = f.ynext(1)[i.template plus<1, direction>()]; } else { s.pp = f[i.template plus<1, direction>()]; } break; case(STAGGER::L2C): if (direction == DIRECTION::YOrthogonal) { - s.mm = f.ydown()[i.template minus<1, direction>()]; + s.mm = f.ynext(-1)[i.template minus<1, direction>()]; } else { s.mm = f[i.template minus<1, direction>()]; } s.m = f[i]; s.c = s.m; if (direction == DIRECTION::YOrthogonal) { - s.p = f.yup()[i.template plus<1, direction>()]; + s.p = f.ynext(1)[i.template plus<1, direction>()]; } else { s.p = f[i.template plus<1, direction>()]; } - if (nGuard == 2) s.pp = f[i.template plus<2, direction>()]; + if (nGuard == 2) { + if (direction == DIRECTION::YOrthogonal) { + s.pp = f.ynext(2)[i.template plus<2, direction>()]; + } else { + s.pp = f[i.template plus<2, direction>()]; + } + } break; } return; diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index ead57190c2..9333032543 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -100,23 +100,8 @@ Field3D::Field3D(const BoutReal val, Mesh* localmesh) : Field(localmesh) { Field3D::~Field3D() { /// Delete the time derivative variable if allocated if (deriv != nullptr) { - // The ddt of the yup/ydown_fields point to the same place as ddt.yup_field - // only delete once - // Also need to check that separate yup_field exists - if ((yup_field != this) && (yup_field != nullptr)) - yup_field->deriv = nullptr; - if ((ydown_field != this) && (ydown_field != nullptr)) - ydown_field->deriv = nullptr; - - // Now delete them as part of the deriv vector delete deriv; } - - if((yup_field != this) && (yup_field != nullptr)) - delete yup_field; - - if((ydown_field != this) && (ydown_field != nullptr)) - delete ydown_field; } void Field3D::allocate() { @@ -146,54 +131,69 @@ Field3D* Field3D::timeDeriv() { void Field3D::splitYupYdown() { TRACE("Field3D::splitYupYdown"); - if((yup_field != this) && (yup_field != nullptr)) +#if CHECK > 2 + if (yup_fields.size() != ydown_fields.size()) { + throw BoutException("Field3D::splitYupYdown: forward/backward parallel slices not in sync.\n" + " This is an internal library error"); + } +#endif + + if (!yup_fields.empty()) { return; + } - // yup_field and ydown_field null - yup_field = new Field3D(fieldmesh); - ydown_field = new Field3D(fieldmesh); + for (int i = 0; i < fieldmesh->ystart; ++i) { + yup_fields.emplace_back(fieldmesh); + ydown_fields.emplace_back(fieldmesh); + } } void Field3D::mergeYupYdown() { TRACE("Field3D::mergeYupYdown"); - - if(yup_field == this && ydown_field == this) - return; - if(yup_field != nullptr){ - delete yup_field; +#if CHECK > 2 + if (yup_fields.size() != ydown_fields.size()) { + throw BoutException("Field3D::mergeYupYdown: forward/backward parallel slices not in sync.\n" + " This is an internal library error"); } +#endif - if(ydown_field != nullptr) { - delete ydown_field; + if (yup_fields.empty() && ydown_fields.empty()) { + return; } - yup_field = this; - ydown_field = this; + yup_fields.clear(); + ydown_fields.clear(); } -Field3D& Field3D::ynext(int dir) { - switch(dir) { - case +1: - return yup(); - case -1: - return ydown(); - default: - throw BoutException("Field3D: Call to ynext with strange direction %d. Only +/-1 currently supported", dir); +const Field3D& Field3D::ynext(int dir) const { +#if CHECK > 0 + // Asked for more than yguards + if (std::abs(dir) > fieldmesh->ystart) { + throw BoutException( + "Field3D: Call to ynext with %d which is more than number of yguards (%d)", dir, + fieldmesh->ystart); } -} +#endif -const Field3D& Field3D::ynext(int dir) const { - switch(dir) { - case +1: - return yup(); - case -1: - return ydown(); - default: - throw BoutException("Field3D: Call to ynext with strange direction %d. Only +/-1 currently supported", dir); + // ynext uses 1-indexing, but yup wants 0-indexing + if (dir > 0) { + return yup(dir - 1); + } else if (dir < 0) { + return ydown(std::abs(dir) - 1); + } else { + return *this; } } +Field3D &Field3D::ynext(int dir) { + // Call the `const` version: need to add `const` to `this` to call + // it, then throw it away after. This is ok because `this` wasn't + // `const` to begin with. + // See Effective C++, Scott Meyers, p23, for a better explanation + return const_cast(static_cast(*this).ynext(dir)); +} + void Field3D::setLocation(CELL_LOC new_location) { AUTO_TRACE(); if (getMesh()->StaggerGrids) { diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index b3ba1cc95f..0a2b8e79ca 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -340,8 +340,8 @@ const Field3D Vpar_Grad_par_LCtoC(const Field3D &v, const Field3D &f, REGION reg result.allocate(); - bool vUseUpDown = (v.hasYupYdown() && ((&v.yup() != &v) || (&v.ydown() != &v))); - bool fUseUpDown = (f.hasYupYdown() && ((&f.yup() != &f) || (&f.ydown() != &f))); + bool vUseUpDown = v.hasYupYdown(); + bool fUseUpDown = f.hasYupYdown(); if (vUseUpDown && fUseUpDown) { // Both v and f have up/down fields diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx index 1d37b1fe7a..d85f7df1e4 100644 --- a/src/mesh/parallel/fci.cxx +++ b/src/mesh/parallel/fci.cxx @@ -42,83 +42,99 @@ #include "parallel_boundary_region.hxx" #include #include -#include // See this for codes +#include #include #include -/** - * Return the sign of val - */ -inline BoutReal sgn(BoutReal val) { return (BoutReal(0) < val) - (val < BoutReal(0)); } +#include -// Calculate all the coefficients needed for the spline interpolation -// dir MUST be either +1 or -1 -FCIMap::FCIMap(Mesh &mesh_in, int dir, bool zperiodic) - : mesh(mesh_in), dir(dir), boundary_mask(mesh_in), corner_boundary_mask(mesh_in), - y_prime(&mesh_in) { +FCIMap::FCIMap(Mesh& mesh, int offset_, BoundaryRegionPar* boundary, bool zperiodic) + : map_mesh(mesh), offset(offset_), boundary_mask(map_mesh), + corner_boundary_mask(map_mesh) { - interp = InterpolationFactory::getInstance()->create(&mesh); - interp->setYOffset(dir); + TRACE("Creating FCIMAP for direction %d", offset); + + if (offset == 0) { + throw BoutException("FCIMap called with offset = 0; You probably didn't mean to do that"); + } + + interp = + std::unique_ptr(InterpolationFactory::getInstance()->create(&map_mesh)); + interp->setYOffset(offset); + + interp_corner = + std::unique_ptr(InterpolationFactory::getInstance()->create(&map_mesh)); + interp_corner->setYOffset(offset); - interp_corner = InterpolationFactory::getInstance()->create(&mesh); - interp_corner->setYOffset(dir); - // Index arrays contain guard cells in order to get subscripts right // x-index of bottom-left grid point - auto i_corner = Tensor(mesh.LocalNx, mesh.LocalNy, mesh.LocalNz); + auto i_corner = Tensor(map_mesh.LocalNx, map_mesh.LocalNy, map_mesh.LocalNz); // z-index of bottom-left grid point - auto k_corner = Tensor(mesh.LocalNx, mesh.LocalNy, mesh.LocalNz); - - Field3D xt_prime(&mesh), zt_prime(&mesh); - Field3D R(&mesh), Z(&mesh); // Real-space coordinates of grid points - Field3D R_prime(&mesh), - Z_prime(&mesh); // Real-space coordinates of forward/backward points - - mesh.get(R, "R", 0.0, false); - mesh.get(Z, "Z", 0.0, false); - - // Load the floating point indices from the grid file - // Future, higher order parallel derivatives could require maps to +/-2 slices - if (dir == +1) { - mesh.get(xt_prime, "forward_xt_prime", 0.0, false); - mesh.get(zt_prime, "forward_zt_prime", 0.0, false); - mesh.get(R_prime, "forward_R", 0.0, false); - mesh.get(Z_prime, "forward_Z", 0.0, false); - boundary = new BoundaryRegionPar("FCI_forward", BNDRY_PAR_FWD, dir, &mesh); - } else if (dir == -1) { - mesh.get(xt_prime, "backward_xt_prime", 0.0, false); - mesh.get(zt_prime, "backward_zt_prime", 0.0, false); - mesh.get(R_prime, "backward_R", 0.0, false); - mesh.get(Z_prime, "backward_Z", 0.0, false); - boundary = new BoundaryRegionPar("FCI_backward", BNDRY_PAR_BKWD, dir, &mesh); - } else { - // Definitely shouldn't be called - throw BoutException("FCIMap called with strange direction: %d. Only +/-1 currently supported.", dir); + auto k_corner = Tensor(map_mesh.LocalNx, map_mesh.LocalNy, map_mesh.LocalNz); + + // Index-space coordinates of forward/backward points + Field3D xt_prime(&map_mesh), zt_prime(&map_mesh); + // Real-space coordinates of grid points + Field3D R(&map_mesh), Z(&map_mesh); + // Real-space coordinates of forward/backward points + Field3D R_prime(&map_mesh), Z_prime(&map_mesh); + + map_mesh.get(R, "R", 0.0, false); + map_mesh.get(Z, "Z", 0.0, false); + + // Get a unique name for a field based on the sign/magnitude of the offset + const auto parallel_slice_field_name = [&](std::string field) -> std::string { + const std::string direction = (offset > 0) ? "forward" : "backward"; + // We only have a suffix for parallel slices beyond the first + // This is for backwards compatibility + const std::string slice_suffix = + (std::abs(offset) > 1) ? "_" + std::to_string(std::abs(offset)) : ""; + return direction + "_" + field + slice_suffix; + }; + + // If we can't read in any of these fields, things will silently not + // work, so best throw + if (map_mesh.get(xt_prime, parallel_slice_field_name("xt_prime"), 0.0, false) != 0) { + throw BoutException("Could not read %s from grid file!\n" + " Either add it to the grid file, or reduce MYG", + parallel_slice_field_name("xt_prime").c_str()); + } + if (map_mesh.get(zt_prime, parallel_slice_field_name("zt_prime"), 0.0, false) != 0) { + throw BoutException("Could not read %s from grid file!\n" + " Either add it to the grid file, or reduce MYG", + parallel_slice_field_name("zt_prime").c_str()); + } + if (map_mesh.get(R_prime, parallel_slice_field_name("R"), 0.0, false) != 0) { + throw BoutException("Could not read %s from grid file!\n" + " Either add it to the grid file, or reduce MYG", + parallel_slice_field_name("R").c_str()); + } + if (map_mesh.get(Z_prime, parallel_slice_field_name("Z"), 0.0, false) != 0) { + throw BoutException("Could not read %s from grid file!\n" + " Either add it to the grid file, or reduce MYG", + parallel_slice_field_name("Z").c_str()); } - // Add the boundary region to the mesh's vector of parallel boundaries - mesh.addBoundaryPar(boundary); - // Cell corners - Field3D xt_prime_corner(&mesh), zt_prime_corner(&mesh); + Field3D xt_prime_corner(&map_mesh), zt_prime_corner(&map_mesh); xt_prime_corner.allocate(); zt_prime_corner.allocate(); - for (int x = mesh.xstart; x <= mesh.xend; x++) { - for (int y = mesh.ystart; y <= mesh.yend; y++) { - for (int z = 0; z < mesh.LocalNz - 1; z++) { + for (int x = map_mesh.xstart; x <= map_mesh.xend; x++) { + for (int y = map_mesh.ystart; y <= map_mesh.yend; y++) { + for (int z = 0; z < map_mesh.LocalNz - 1; z++) { // Point interpolated from (x+1/2, z+1/2) if ((xt_prime(x, y, z) < 0.0) || (xt_prime(x + 1, y, z) < 0.0) || (xt_prime(x + 1, y, z + 1) < 0.0) || (xt_prime(x, y, z + 1) < 0.0)) { // Hit a boundary corner_boundary_mask(x, y, z) = true; - + xt_prime_corner(x, y, z) = -1.0; zt_prime_corner(x, y, z) = -1.0; continue; } - + xt_prime_corner(x, y, z) = 0.25 * (xt_prime(x, y, z) + xt_prime(x + 1, y, z) + xt_prime(x, y, z + 1) + xt_prime(x + 1, y, z + 1)); @@ -129,19 +145,19 @@ FCIMap::FCIMap(Mesh &mesh_in, int dir, bool zperiodic) } } } - + interp_corner->setMask(corner_boundary_mask); interp_corner->calcWeights(xt_prime_corner, zt_prime_corner); - + interp->calcWeights(xt_prime, zt_prime); - - int ncz = mesh.LocalNz; + + int ncz = map_mesh.LocalNz; BoutReal t_x, t_z; - Coordinates &coord = *(mesh.getCoordinates()); + Coordinates &coord = *(map_mesh.getCoordinates()); - for (int x = mesh.xstart; x <= mesh.xend; x++) { - for (int y = mesh.ystart; y <= mesh.yend; y++) { + for (int x = map_mesh.xstart; x <= map_mesh.xend; x++) { + for (int y = map_mesh.ystart; y <= map_mesh.yend; y++) { for (int z = 0; z < ncz; z++) { // The integer part of xt_prime, zt_prime are the indices of the cell @@ -198,7 +214,7 @@ FCIMap::FCIMap(Mesh &mesh_in, int dir, bool zperiodic) dR_dz = R(x, y, z + 1) - R(x, y, z); dZ_dz = Z(x, y, z + 1) - Z(x, y, z); - } else if (z == mesh.LocalNz - 1) { + } else if (z == map_mesh.LocalNz - 1) { dR_dz = R(x, y, z) - R(x, y, z - 1); dZ_dz = Z(x, y, z) - Z(x, y, z - 1); @@ -215,8 +231,8 @@ FCIMap::FCIMap(Mesh &mesh_in, int dir, bool zperiodic) // Invert 2x2 matrix to get change in index BoutReal dx = (dZ_dz * dR - dR_dz * dZ) / det; BoutReal dz = (dR_dx * dZ - dZ_dx * dR) / det; - boundary->add_point(x, y, z, - x + dx, y + 0.5*dir, z + dz, // Intersection point in local index space + boundary->add_point(x, y, z, + x + dx, y + 0.5*offset, z + dz, // Intersection point in local index space 0.5*coord.dy(x,y), //sqrt( SQ(dR) + SQ(dZ) ), // Distance to intersection PI // Right-angle intersection ); @@ -237,12 +253,14 @@ FCIMap::FCIMap(Mesh &mesh_in, int dir, bool zperiodic) interp->setMask(boundary_mask); } -const Field3D FCIMap::integrate(Field3D &f) const { +Field3D FCIMap::integrate(Field3D &f) const { TRACE("FCIMap::integrate"); - + + ASSERT3(&map_mesh == f.getMesh()); + // Cell centre values Field3D centre = interp->interpolate(f); - + // Cell corner values (x+1/2, z+1/2) Field3D corner = interp_corner->interpolate(f); @@ -250,27 +268,27 @@ const Field3D FCIMap::integrate(Field3D &f) const { result.allocate(); result.setLocation(f.getLocation()); - int nz = mesh.LocalNz; - - for(int x = mesh.xstart; x <= mesh.xend; x++) { - for(int y = mesh.ystart; y <= mesh.yend; y++) { - - int ynext = y+dir; - + int nz = map_mesh.LocalNz; + + for(int x = map_mesh.xstart; x <= map_mesh.xend; x++) { + for(int y = map_mesh.ystart; y <= map_mesh.yend; y++) { + + int ynext = y+offset; + for(int z = 0; z < nz; z++) { if (boundary_mask(x,y,z)) continue; - + int zm = z - 1; if (z == 0) { zm = nz-1; } - + BoutReal f_c = centre(x,ynext,z); - + if (corner_boundary_mask(x, y, z) || corner_boundary_mask(x - 1, y, z) || corner_boundary_mask(x, y, zm) || corner_boundary_mask(x - 1, y, zm) || - (x == mesh.xstart)) { + (x == map_mesh.xstart)) { // One of the corners leaves the domain. // Use the cell centre value, since boundary conditions are not // currently applied to corners. @@ -295,24 +313,26 @@ const Field3D FCIMap::integrate(Field3D &f) const { return result; } -void FCITransform::calcYUpDown(Field3D &f) { +void FCITransform::calcYUpDown(Field3D& f) { TRACE("FCITransform::calcYUpDown"); // Ensure that yup and ydown are different fields f.splitYupYdown(); // Interpolate f onto yup and ydown fields - f.ynext(forward_map.dir) = forward_map.interpolate(f); - f.ynext(backward_map.dir) = backward_map.interpolate(f); + for (const auto& map : field_line_maps) { + f.ynext(map.offset) = map.interpolate(f); + } } -void FCITransform::integrateYUpDown(Field3D &f) { +void FCITransform::integrateYUpDown(Field3D& f) { TRACE("FCITransform::integrateYUpDown"); - + // Ensure that yup and ydown are different fields f.splitYupYdown(); // Integrate f onto yup and ydown fields - f.ynext(forward_map.dir) = forward_map.integrate(f); - f.ynext(backward_map.dir) = backward_map.integrate(f); + for (const auto& map : field_line_maps) { + f.ynext(map.offset) = map.integrate(f); + } } diff --git a/src/mesh/parallel/fci.hxx b/src/mesh/parallel/fci.hxx index d9cfe5af0c..37a87a96e6 100644 --- a/src/mesh/parallel/fci.hxx +++ b/src/mesh/parallel/fci.hxx @@ -32,44 +32,59 @@ #include #include -/*! - * Field line map - contains the coefficients for interpolation - */ -class FCIMap { - /// Interpolation object - Interpolation *interp; // Cell centre - Interpolation *interp_corner; // Cell corner at (x+1, z+1) +#include +#include - Mesh& mesh; - /// Private constructor - must be initialised with mesh - FCIMap(); -public: - /// dir MUST be either +1 or -1 - FCIMap(Mesh& mesh_in, int dir, bool zperiodic); +/// Field line map - contains the coefficients for interpolation +class FCIMap { + /// Interpolation objects + std::unique_ptr interp; // Cell centre + std::unique_ptr interp_corner; // Cell corner at (x+1, z+1) - int dir; /**< Direction of map */ +public: + FCIMap() = delete; + FCIMap(Mesh& mesh, int offset, BoundaryRegionPar* boundary, bool zperiodic); - BoutMask boundary_mask; /**< boundary mask - has the field line left the domain */ - BoutMask corner_boundary_mask; ///< If any of the integration area has left the domain - - Field3D y_prime; /**< distance to intersection with boundary */ + // The mesh this map was created on + Mesh& map_mesh; - BoundaryRegionPar* boundary; /**< boundary region */ + /// Direction of map + const int offset; - const Field3D interpolate(Field3D &f) const { return interp->interpolate(f); } + /// boundary mask - has the field line left the domain + BoutMask boundary_mask; + /// If any of the integration area has left the domain + BoutMask corner_boundary_mask; + + Field3D interpolate(Field3D& f) const { + ASSERT3(&map_mesh == f.getMesh()); + return interp->interpolate(f); + } - const Field3D integrate(Field3D &f) const; + Field3D integrate(Field3D &f) const; }; -/*! - * Flux Coordinate Independent method for parallel derivatives - */ + +/// Flux Coordinate Independent method for parallel derivatives class FCITransform : public ParallelTransform { public: - FCITransform(Mesh &mesh, bool zperiodic = true) - : mesh(mesh), forward_map(mesh, +1, zperiodic), backward_map(mesh, -1, zperiodic), - zperiodic(zperiodic) {} + FCITransform() = delete; + FCITransform(Mesh& mesh, bool zperiodic = true) { + + auto forward_boundary = new BoundaryRegionPar("FCI_forward", BNDRY_PAR_FWD, +1, &mesh); + auto backward_boundary = new BoundaryRegionPar("FCI_backward", BNDRY_PAR_BKWD, -1, &mesh); + + // Add the boundary region to the mesh's vector of parallel boundaries + mesh.addBoundaryPar(forward_boundary); + mesh.addBoundaryPar(backward_boundary); + + field_line_maps.reserve(mesh.ystart * 2); + for (int offset = 1; offset < mesh.ystart + 1; ++offset) { + field_line_maps.emplace_back(mesh, offset, forward_boundary, zperiodic); + field_line_maps.emplace_back(mesh, -offset, backward_boundary, zperiodic); + } + } void calcYUpDown(Field3D &f) override; @@ -83,18 +98,11 @@ public: throw BoutException("FCI method cannot transform into field aligned grid"); } - bool canToFromFieldAligned() override{ - return false; - } -private: - FCITransform(); - - Mesh& mesh; + bool canToFromFieldAligned() override { return false; } - FCIMap forward_map; /**< FCI map for field lines in +ve y */ - FCIMap backward_map; /**< FCI map for field lines in -ve y */ - - bool zperiodic; /**< Is the z-direction periodic? */ +private: + /// FCI maps for each of the parallel slices + std::vector field_line_maps; }; #endif // __FCITRANSFORM_H__ diff --git a/src/mesh/parallel/identity.cxx b/src/mesh/parallel/identity.cxx new file mode 100644 index 0000000000..437d04178b --- /dev/null +++ b/src/mesh/parallel/identity.cxx @@ -0,0 +1,11 @@ +#include "bout/paralleltransform.hxx" +#include "bout/mesh.hxx" + +void ParallelTransformIdentity::calcYUpDown(Field3D& f) { + f.splitYupYdown(); + + for (int i = 0; i < f.getMesh()->ystart; ++i) { + f.yup(i) = f; + f.ydown(i) = f; + } +} diff --git a/src/mesh/parallel/makefile b/src/mesh/parallel/makefile index 9c4817b008..2e522318de 100644 --- a/src/mesh/parallel/makefile +++ b/src/mesh/parallel/makefile @@ -2,7 +2,7 @@ BOUT_TOP = ../../.. DIRS = -SOURCEC = shiftedmetric.cxx fci.cxx +SOURCEC = shiftedmetric.cxx fci.cxx identity.cxx TARGET = lib include $(BOUT_TOP)/make.config diff --git a/src/mesh/parallel/shiftedmetric.cxx b/src/mesh/parallel/shiftedmetric.cxx index 55eee670bd..d6171fc117 100644 --- a/src/mesh/parallel/shiftedmetric.cxx +++ b/src/mesh/parallel/shiftedmetric.cxx @@ -58,21 +58,13 @@ void ShiftedMetric::cachePhases() { fromAlignedPhs.resize(mesh.LocalNx); toAlignedPhs.resize(mesh.LocalNx); - yupPhs.resize(mesh.LocalNx); - ydownPhs.resize(mesh.LocalNx); - for (int jx = 0; jx < mesh.LocalNx; jx++) { fromAlignedPhs[jx].resize(mesh.LocalNy); toAlignedPhs[jx].resize(mesh.LocalNy); - yupPhs[jx].resize(mesh.LocalNy); - ydownPhs[jx].resize(mesh.LocalNy); for (int jy = 0; jy < mesh.LocalNy; jy++) { fromAlignedPhs[jx][jy].resize(nmodes); toAlignedPhs[jx][jy].resize(nmodes); - - yupPhs[jx][jy].resize(nmodes); - ydownPhs[jx][jy].resize(nmodes); } } @@ -89,48 +81,53 @@ void ShiftedMetric::cachePhases() { } } - // Yup/Ydown phases -- note we don't shift in the boundaries/guards - for (int jx = 0; jx < mesh.LocalNx; jx++) { - for (int jy = mesh.ystart; jy <= mesh.yend; jy++) { - BoutReal yupShift = zShift(jx, jy) - zShift(jx, jy + 1); - BoutReal ydownShift = zShift(jx, jy) - zShift(jx, jy - 1); - - for (int jz = 0; jz < nmodes; jz++) { - BoutReal kwave = jz * 2.0 * PI / zlength; // wave number is 1/[rad] - - yupPhs[jx][jy][jz] = dcomplex(cos(kwave * yupShift), -sin(kwave * yupShift)); - ydownPhs[jx][jy][jz] = - dcomplex(cos(kwave * ydownShift), -sin(kwave * ydownShift)); - } - } + // Allocate space for parallel slice caches: y-guard cells in each + // direction + parallel_slice_phases.resize(mesh.ystart * 2); + + // Careful with the indices/offsets! Offsets are 1-indexed (as 0 + // would be the original slice), and Mesh::ystart is the number of + // guard cells. The parallel slice vector stores the offsets as + // {+1, ..., +n, -1, ..., -n} + // Once parallel_slice_phases is initialised though, each element + // stores its phase and offset, so we don't need to faff about after + // this + for (int i = 0; i < mesh.ystart; ++i) { + // NOTE: std::vector constructor here takes a **copy** of the + // Array! We *must* call `Array::ensureUnique` on each element + // before using it! + parallel_slice_phases[i].phase_shift = + arr3Dvec(mesh.LocalNx, + std::vector>(mesh.LocalNy, Array(nmodes))); + parallel_slice_phases[i].y_offset = i + 1; + + // Backwards parallel slices + parallel_slice_phases[mesh.ystart + i].phase_shift = + arr3Dvec(mesh.LocalNx, + std::vector>(mesh.LocalNy, Array(nmodes))); + parallel_slice_phases[mesh.ystart + i].y_offset = -(i + 1); } -} -/*! - * Calculate the Y up and down fields - */ -void ShiftedMetric::calcYUpDown(Field3D &f) { - f.splitYupYdown(); - - Field3D& yup = f.yup(); - yup.allocate(); + // Parallel slice phases -- note we don't shift in the boundaries/guards + for (auto& slice : parallel_slice_phases) { + for (int jx = 0; jx < mesh.LocalNx; jx++) { + for (int jy = mesh.ystart; jy <= mesh.yend; jy++) { - for(int jx=0;jx& phs, +void ShiftedMetric::shiftZ(const BoutReal* in, const Array& phs, BoutReal* out) const { int nmodes = mesh.LocalNz / 2 + 1; @@ -187,6 +184,67 @@ void ShiftedMetric::shiftZ(const BoutReal* in, const std::vector& phs, irfft(&cmplx[0], mesh.LocalNz, out); // Reverse FFT } + +void ShiftedMetric::calcYUpDown(Field3D& f) { + + auto results = shiftZ(f, parallel_slice_phases); + + ASSERT3(results.size() == parallel_slice_phases.size()); + + f.splitYupYdown(); + + for (std::size_t i = 0; i < results.size(); ++i) { + f.ynext(parallel_slice_phases[i].y_offset) = std::move(results[i]); + } +} + +std::vector +ShiftedMetric::shiftZ(const Field3D& f, + const std::vector& phases) const { + + const int nmodes = mesh.LocalNz / 2 + 1; + + // FFT in Z of input field at each (x, y) point + arr3Dvec f_fft(mesh.LocalNx, + std::vector>(mesh.LocalNy, Array(nmodes))); + + for (int jx = 0; jx < mesh.LocalNx; jx++) { + for (int jy = 0; jy < mesh.LocalNy; jy++) { + f_fft[jx][jy].ensureUnique(); + rfft(f(jx, jy), mesh.LocalNz, f_fft[jx][jy].begin()); + } + } + + std::vector results{}; + + for (auto& phase : phases) { + // In C++17 std::vector::emplace_back returns a reference, which + // would be very useful here! + results.emplace_back(&mesh); + auto& current_result = results.back(); + current_result.allocate(); + current_result.setLocation(f.getLocation()); + + for (int jx = 0; jx < mesh.LocalNx; jx++) { + for (int jy = mesh.ystart; jy <= mesh.yend; jy++) { + + // Deep copy the FFT'd field + Array shifted_temp(f_fft[jx][jy + phase.y_offset]); + shifted_temp.ensureUnique(); + + for (int jz = 1; jz < nmodes; ++jz) { + shifted_temp[jz] *= phase.phase_shift[jx][jy][jz]; + } + + irfft(shifted_temp.begin(), mesh.LocalNz, + current_result(jx, jy + phase.y_offset)); + } + } + } + + return results; +} + //Old approach retained so we can still specify a general zShift const Field3D ShiftedMetric::shiftZ(const Field3D &f, const Field2D &zangle) const { ASSERT1(&mesh == f.getMesh()); diff --git a/tests/MMS/.gitignore b/tests/MMS/.gitignore index 2946406b9b..c8f94147b3 100644 --- a/tests/MMS/.gitignore +++ b/tests/MMS/.gitignore @@ -30,6 +30,9 @@ BOUT.settings /spatial/d2dx2/test_d2dx2 /spatial/d2dz2/test_d2dz2 /spatial/diffusion/diffusion +/spatial/fci/fci_mms +/spatial/fci/fci.grid.nc +/spatial/fci/fci_mms.pkl /time/time /tokamak/tokamak /tokamak/tokamak.pkl diff --git a/tests/MMS/spatial/fci/data/BOUT.inp b/tests/MMS/spatial/fci/data/BOUT.inp new file mode 100644 index 0000000000..cf2208e3b3 --- /dev/null +++ b/tests/MMS/spatial/fci/data/BOUT.inp @@ -0,0 +1,23 @@ +grid = fci.grid.nc + +input = sin(y - 2*z) + sin(y - z) + +solution = (6.28318530717959*(0.01*x + 0.045)*(-2*cos(y - 2*z) - cos(y - z)) + 0.628318530717959*cos(y - 2*z) + 0.628318530717959*cos(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) + +MXG = 1 +NXPE = 1 + +[mesh] +paralleltransform = fci +symmetricglobalx = true + +[mesh:ddy] +first = C2 +second = C2 + +[fci] +y_periodic = true +z_periodic = true + +[interpolation] +type = lagrange4pt diff --git a/tests/MMS/spatial/fci/fci_mms.cxx b/tests/MMS/spatial/fci/fci_mms.cxx new file mode 100644 index 0000000000..dc250d55b7 --- /dev/null +++ b/tests/MMS/spatial/fci/fci_mms.cxx @@ -0,0 +1,28 @@ +#include "bout.hxx" +#include "derivs.hxx" +#include "field_factory.hxx" + +int main(int argc, char** argv) { + BoutInitialise(argc, argv); + + Field3D input{FieldFactory::get()->create3D("input", Options::getRoot(), mesh)}; + Field3D solution{FieldFactory::get()->create3D("solution", Options::getRoot(), mesh)}; + + // Communicate to calculate parallel transform + mesh->communicate(input); + + Field3D result{Grad_par(input)}; + Field3D error{result - solution}; + BoutReal l_2{sqrt(mean(SQ(error), true, RGN_NOBNDRY))}; + BoutReal l_inf{max(abs(error), true, RGN_NOBNDRY)}; + + SAVE_ONCE6(input, solution, result, error, l_2, l_inf); + + for (int slice = 1; slice < mesh->ystart; ++slice) { + SAVE_ONCE2(input.ynext(-slice), input.ynext(slice)); + } + + dump.write(); + + BoutFinalise(); +} diff --git a/tests/MMS/spatial/fci/makefile b/tests/MMS/spatial/fci/makefile new file mode 100644 index 0000000000..88ba6c77e7 --- /dev/null +++ b/tests/MMS/spatial/fci/makefile @@ -0,0 +1,6 @@ + +BOUT_TOP = ../../../.. + +SOURCEC = fci_mms.cxx + +include $(BOUT_TOP)/make.config diff --git a/tests/integrated/test-fci-slab/mms.py b/tests/MMS/spatial/fci/mms.py old mode 100644 new mode 100755 similarity index 50% rename from tests/integrated/test-fci-slab/mms.py rename to tests/MMS/spatial/fci/mms.py index df00c5d68d..477c605e84 --- a/tests/integrated/test-fci-slab/mms.py +++ b/tests/MMS/spatial/fci/mms.py @@ -1,3 +1,4 @@ +#!/usr/bin/env python3 # # Generate manufactured solution and sources for FCI test # @@ -11,9 +12,7 @@ from math import pi -f = sin(y - z) + cos(t)*sin(y - 2*z) - -g = cos(y - z) - cos(t)*sin(y - 2*z) +f = sin(y - z) + sin(y - 2*z) Lx = 0.1 Ly = 10. @@ -26,22 +25,11 @@ Bpx = Bp + (x-0.5)*Lx * Bpprime # Note: x in range [0,1] B = sqrt(Bpx**2 + Bt**2) -def FCI_Grad_par(f): +def FCI_ddy(f): return ( Bt * diff(f, y)*2.*pi/Ly + Bpx * diff(f, z)*2.*pi/Lz ) / B ############################################ # Equations solved -dfdt = FCI_Grad_par(g) -dgdt = FCI_Grad_par(f) - -# Loop over variables and print solution, source etc. -for v, dvdt, name in [ (f, dfdt, "f"), (g, dgdt, "g") ]: - # Calculate source - S = diff(v, t) - dvdt - - print("\n["+name+"]") - print("solution = "+exprToStr(v)) - print("\nsource = "+exprToStr(S)) - print("\nbndry_par_all = parallel_dirichlet("+name+":solution)") - +print("input = " + exprToStr(f)) +print("solution = " + exprToStr(FCI_ddy(f))) diff --git a/tests/MMS/spatial/fci/runtest b/tests/MMS/spatial/fci/runtest new file mode 100755 index 0000000000..99ec92edc2 --- /dev/null +++ b/tests/MMS/spatial/fci/runtest @@ -0,0 +1,176 @@ +#!/usr/bin/env python3 +# +# Python script to run and analyse MMS test +# +from __future__ import division +from __future__ import print_function + +from boututils.run_wrapper import shell_safe, launch_safe, getmpirun +from boutdata.collect import collect + +from numpy import array, log, polyfit, linspace, arange + +import pickle + +from sys import stdout + +import zoidberg as zb + +nx = 3 # Not changed for these tests + +# Resolution in y and z +nlist = [8, 16, 32, 64, 128] + +# Number of parallel slices (in each direction) +nslices = [1, 2] + +directory = "data" + +nproc = 2 +mthread = 2 + +MPIRUN = getmpirun() + +success = True + +error_2 = {} +error_inf = {} +method_orders = {} + +# Run with periodic Y? +yperiodic = True + +failures = [] + +print("Making fci MMS test") +shell_safe("make > make.log") + +for nslice in nslices: + error_2[nslice] = [] + error_inf[nslice] = [] + + # Which central difference scheme to use and its expected order + order = nslice * 2 + method_orders[nslice] = { + "name": "C{}".format(order), + "order": order + } + + for n in nlist: + # Define the magnetic field using new poloidal gridding method + # Note that the Bz and Bzprime parameters here must be the same as in mms.py + field = zb.field.Slab(Bz=0.05, Bzprime=0.1) + # Create rectangular poloidal grids + poloidal_grid = zb.poloidal_grid.RectangularPoloidalGrid(nx, n, 0.1, 1.) + # Set the ylength and y locations + ylength = 10. + + if yperiodic: + ycoords = linspace(0.0, ylength, n, endpoint=False) + else: + # Doesn't include the end points + ycoords = (arange(n) + 0.5)*ylength/float(n) + + # Create the grid + grid = zb.grid.Grid(poloidal_grid, ycoords, ylength, yperiodic=yperiodic) + # Make and write maps + maps = zb.make_maps(grid, field, nslice=nslice, quiet=True) + zb.write_maps(grid, field, maps, new_names=False, metric2d=True, quiet=True) + + args = (" MZ={} MYG={} fci:y_periodic={} mesh:ddy:first={}" + .format(n, nslice, yperiodic, method_orders[nslice]["name"])) + + # Command to run + cmd = "./fci_mms "+args + + print("Running command: "+cmd) + + # Launch using MPI + s, out = launch_safe(cmd, runcmd=MPIRUN, nproc=nproc, mthread=mthread, pipe=True) + + # Save output to log file + with open("run.log."+str(n), "w") as f: + f.write(out) + + if s: + print("Run failed!\nOutput was:\n") + print(out) + exit(s) + + # Collect data + l_2 = collect("l_2", tind=[1, 1], info=False, + path=directory, xguards=False, yguards=False) + l_inf = collect("l_inf", tind=[1, 1], info=False, + path=directory, xguards=False, yguards=False) + + error_2[nslice].append(l_2) + error_inf[nslice].append(l_inf) + + print("Errors : l-2 {:f} l-inf {:f}".format(l_2, l_inf)) + + dx = 1. / array(nlist) + + # Calculate convergence order + fit = polyfit(log(dx), log(error_2[nslice]), 1) + order = fit[0] + stdout.write("Convergence order = {:f} (fit)".format(order)) + + order = log(error_2[nslice][-2]/error_2[nslice][-1])/log(dx[-2]/dx[-1]) + stdout.write(", {:f} (small spacing)".format(order)) + + # Should be close to the expected order + if order > method_orders[nslice]["order"] * 0.95: + print("............ PASS\n") + else: + print("............ FAIL\n") + success = False + failures.append(method_orders[nslice]["name"]) + + +with open("fci_mms.pkl", "wb") as output: + pickle.dump(nlist, output) + for nslice in nslices: + pickle.dump(error_2[nslice], output) + pickle.dump(error_inf[nslice], output) + +# Do we want to show the plot as well as save it to file. +showPlot = True + +if False: + try: + # Plot using matplotlib if available + import matplotlib.pyplot as plt + + fig, ax = plt.subplots(1, 1) + + for nslice in nslices: + ax.plot(dx, error_2[nslice], '-', + label="{} $l_2$".format(method_orders[nslice]["name"])) + ax.plot(dx, error_inf[nslice], '--', + label="{} $l_\inf$".format(method_orders[nslice]["name"])) + ax.legend(loc="upper left") + ax.grid() + ax.set_yscale('log') + ax.set_xscale('log') + ax.set_title('error scaling') + ax.set_xlabel(r'Mesh spacing $\delta x$') + ax.set_ylabel("Error norm") + + plt.savefig("fci_mms.pdf") + + print("Plot saved to fci_mms.pdf") + + if showPlot: + plt.show() + plt.close() + except ImportError: + print("No matplotlib") + +if success: + print("All tests passed") + exit(0) +else: + print("Some tests failed:") + for failure in failures: + print("\t" + failure) + exit(1) diff --git a/tests/integrated/test-fci-slab/data/BOUT.inp b/tests/integrated/test-fci-slab/data/BOUT.inp deleted file mode 100644 index 218946e676..0000000000 --- a/tests/integrated/test-fci-slab/data/BOUT.inp +++ /dev/null @@ -1,26 +0,0 @@ -grid = fci.grid.nc -#grid = simple_test.nc - -nout = 50 -timestep = 0.2 - -MZ = 64 - -[mesh] -paralleltransform = fci - -symmetricglobalx = true - -[interpolation] -type=lagrange4pt - -[fci] -y_periodic = false -z_periodic = false - -[f] -scale = 1.0 -function = cos(y-z) - -bndry_par_yup = parallel_dirichlet(0.0) -bndry_par_ydown = parallel_dirichlet(0.0) diff --git a/tests/integrated/test-fci-slab/fci.py b/tests/integrated/test-fci-slab/fci.py deleted file mode 100644 index 1c505ffa20..0000000000 --- a/tests/integrated/test-fci-slab/fci.py +++ /dev/null @@ -1,110 +0,0 @@ -from __future__ import division -from builtins import object -from past.utils import old_div -import numpy as np -from math import pi -from scipy.integrate import odeint -import boututils.datafile as bdata -from boutdata.input import transform3D - -# Parameters -nx = 34 -########## y is toroidal! -ny = 64 -########## z is poloidal! -nz = 64 - -Lx = 0.1 # Radial domain size [m] -Ltor = 10. # "Toroidal" length [m] -Lpol = 1. # "Poloidal" length [m] - -delta_x = old_div(Lx,(nx)) -delta_pol = old_div(Lpol,(nz)) -delta_tor = old_div(Ltor,(ny)) - -Bt = 1.0 # Magnetic field [T] -Bp = 0.1 # Poloidal field at the middle of the domain [T] -Bpprime = 1.0 # Bp gradient [T/m] Bp(x) = Bp + Bpprime * x - -# Coord arrays -x = np.linspace(0,Lx,nx) -y = np.linspace(0,Ltor,ny) -z = np.linspace(0,Lpol,nz,endpoint=False) - -############################################################ - -# Effective major radius -R = old_div(Ltor, (2.*pi)) - -# Set poloidal magnetic field - -Bpx = Bp + (x-old_div(Lx,2)) * Bpprime - -Bpxy = np.transpose(np.resize(Bpx, (nz, ny, nx)), (2,1,0)) - -Bxy = np.sqrt(Bpxy**2 + Bt**2) - -############################################################ - -class Mappoint(object): - def __init__(self, xt, zt): - self.xt = xt - self.zt = zt - - self.xt_prime = old_div(xt,delta_x) - self.zt_prime = old_div(zt,delta_pol) - -def unroll_map_coeff(map_list, coeff): - coeff_array = np.transpose(np.resize(np.array([getattr(f, coeff) for f in map_list]).reshape( (nx,nz) ), (ny, nx, nz) ), (1, 0, 2) ) - return coeff_array - -def b_field(vector, y): - x0 = 0.05 # Centre of box, where bz = 0. - x, z = vector; - bx = 0. - bz = Bp + (x-x0) * Bpprime - - return [bx, bz] - -def field_line_tracer(direction, map_list): - - result = np.zeros( (nx, nz, 2) ) - - for i in np.arange(0,nx): - for k in np.arange(0,nz): - result[i,k,:] = odeint(b_field, [x[i], z[k]], [0, delta_tor*direction])[1,:] - result[i,k,1] = np.mod(result[i,k,1], Lpol) - - map_list.append(Mappoint(result[i,k,0],result[i,k,1])) - - return result - -if __name__ == "__main__": - - forward_map = [] - forward_coords = field_line_tracer(+1, forward_map) - backward_map = [] - backward_coords = field_line_tracer(-1, backward_map) - - X,Y = np.meshgrid(x,y,indexing='ij') - x0 = 0.5 - g_22 = np.sqrt(((Bp + (X-x0) * Lx * Bpprime)**2 + 1)) - - with bdata.DataFile('fci.grid.nc', write=True, create=True) as f: - f.write('nx', nx) - f.write('ny', ny) - f.write('nz', nz) - f.write("dx", delta_x) - f.write("dy", delta_tor) - f.write("g_22", g_22) - f.write("Bxy", transform3D(Bxy)) - - xt_prime = unroll_map_coeff(forward_map, 'xt_prime') - f.write('forward_xt_prime', transform3D(xt_prime)) - zt_prime = unroll_map_coeff(forward_map, 'zt_prime') - f.write('forward_zt_prime', transform3D(zt_prime)) - - xt_prime = unroll_map_coeff(backward_map, 'xt_prime') - f.write('backward_xt_prime', transform3D(xt_prime)) - zt_prime = unroll_map_coeff(backward_map, 'zt_prime') - f.write('backward_zt_prime', transform3D(zt_prime)) diff --git a/tests/integrated/test-fci-slab/fci_slab.cxx b/tests/integrated/test-fci-slab/fci_slab.cxx deleted file mode 100644 index 60a6bfd4ca..0000000000 --- a/tests/integrated/test-fci-slab/fci_slab.cxx +++ /dev/null @@ -1,53 +0,0 @@ -#include -#include -#include - -class FCISlab : public PhysicsModel { -public: - - // We need to initialise the FCI object with the mesh - FCISlab() {} - - int init(bool UNUSED(restarting)) { - - D = 10; - - Coordinates *coord = mesh->getCoordinates(); - - mesh->get(coord->g_22, "g_22"); - - coord->geometry(); - - solver->add(f, "f"); - solver->add(g, "g"); - - f.applyBoundary("dirichlet"); - g.applyBoundary("dirichlet"); - - return 0; - } - - int rhs(BoutReal time); - -private: - Field3D f, g; - - BoutReal D; -}; - -BOUTMAIN(FCISlab); - -int FCISlab::rhs(BoutReal time) { - mesh->communicate(f,g); - - Coordinates *coord = mesh->getCoordinates(); - - f.applyParallelBoundary(time); - g.applyParallelBoundary(time); - - ddt(f) = Grad_par(g) + D*SQ(coord->dy)*Grad2_par2(f); - - ddt(g) = Grad_par(f) + D*SQ(coord->dy)*Grad2_par2(g); - - return 0; -} diff --git a/tests/integrated/test-fci-slab/generate.py b/tests/integrated/test-fci-slab/generate.py deleted file mode 100644 index 03767db3a6..0000000000 --- a/tests/integrated/test-fci-slab/generate.py +++ /dev/null @@ -1,126 +0,0 @@ -from __future__ import division -from builtins import object -from past.utils import old_div -# -# Routines to generate slab meshes for FCI -# - -import numpy as np -from math import pi -from scipy.integrate import odeint -import boututils.datafile as bdata - - -def slab(nx, ny, nz, - filename="fci.grid.nc", - Lx=0.1, Ly=10., Lz = 1., - Bt=1.0, Bp = 0.1, Bpprime = 1.0): - """ - nx - Number of radial points - ny - Number of toroidal points (NOTE: Different to BOUT++ standard) - nz - Number of poloidal points - - Lx - Radial domain size [m] - Ly - Toroidal domain size [m] - Lz - Poloidal domain size [m] - - Bt - Toroidal magnetic field [T] - Bp - Poloidal magnetic field [T] - Bpprime - Gradient of Bp [T/m] Bp(x) = Bp + Bpprime * x - """ - - MXG = 2 - - # Make sure input types are sane - nx = int(nx) - ny = int(ny) - nz = int(nz) - - Lx = float(Lx) - Ly = float(Ly) - Lz = float(Lz) - - delta_x = old_div(Lx,(nx-2.*MXG)) - delta_pol = old_div(Lz,(nz)) - delta_tor = old_div(Ly,(ny)) - - # Coord arrays - x = Lx * (np.arange(nx) - MXG + 0.5)/(nx - 2.*MXG) # 0 and 1 half-way between cells - y = np.linspace(0,Ly,ny) - z = np.linspace(0,Lz,nz,endpoint=False) - - ############################################################ - - # Effective major radius - R = old_div(Ly, (2.*pi)) - - # Set poloidal magnetic field - - Bpx = Bp + (x-old_div(Lx,2)) * Bpprime - - Bpxy = np.transpose(np.resize(Bpx, (nz, ny, nx)), (2,1,0)) - - Bxy = np.sqrt(Bpxy**2 + Bt**2)[:,:,0] - - class Mappoint(object): - def __init__(self, xt, zt): - self.xt = xt - self.zt = zt - - self.xt_prime = old_div(xt,delta_x) + MXG - 0.5 - self.zt_prime = old_div(zt,delta_pol) - - def unroll_map_coeff(map_list, coeff): - coeff_array = np.transpose(np.resize(np.array([getattr(f, coeff) for f in map_list]).reshape( (nx,nz) ), (ny, nx, nz) ), (1, 0, 2) ) - return coeff_array - - def b_field(vector, y): - x0 = old_div(Lx,2.) # Centre of box, where bz = 0. - x, z = vector; - bx = 0. - bz = Bp + (x-x0) * Bpprime - - return [bx, bz] - - def field_line_tracer(direction, map_list): - - result = np.zeros( (nx, nz, 2) ) - - for i in np.arange(0,nx): - for k in np.arange(0,nz): - result[i,k,:] = odeint(b_field, [x[i], z[k]], [0, delta_tor*direction])[1,:] - map_list.append(Mappoint(result[i,k,0],result[i,k,1])) - - return result - - forward_map = [] - forward_coords = field_line_tracer(+1, forward_map) - backward_map = [] - backward_coords = field_line_tracer(-1, backward_map) - - X,Y = np.meshgrid(x,y,indexing='ij') - x0 = 0.5 - g_22 = old_div(((Bp + (X-x0) * Lx * Bpprime)**2 + Bt**2), Bt**2) - - with bdata.DataFile(filename, write=True, create=True) as f: - f.write('nx', nx) - f.write('ny', ny) - f.write('nz', nz) - f.write("dx", delta_x) - f.write("dy", delta_tor) - f.write("g_22", g_22) - f.write("Bxy", (Bxy)) - - xt_prime = unroll_map_coeff(forward_map, 'xt_prime') - f.write('forward_xt_prime', (xt_prime)) - zt_prime = unroll_map_coeff(forward_map, 'zt_prime') - f.write('forward_zt_prime', (zt_prime)) - - xt_prime = unroll_map_coeff(backward_map, 'xt_prime') - f.write('backward_xt_prime', (xt_prime)) - zt_prime = unroll_map_coeff(backward_map, 'zt_prime') - f.write('backward_zt_prime', (zt_prime)) - - -if __name__ == "__main__": - slab(34, 64, 64, filename="fci.grid.nc") diff --git a/tests/integrated/test-fci-slab/makefile b/tests/integrated/test-fci-slab/makefile deleted file mode 100644 index b0fbe93385..0000000000 --- a/tests/integrated/test-fci-slab/makefile +++ /dev/null @@ -1,6 +0,0 @@ - -BOUT_TOP = ../../.. - -SOURCEC = fci_slab.cxx - -include $(BOUT_TOP)/make.config diff --git a/tests/integrated/test-fci-slab/mms/BOUT.inp b/tests/integrated/test-fci-slab/mms/BOUT.inp deleted file mode 100644 index 3c1e5f9559..0000000000 --- a/tests/integrated/test-fci-slab/mms/BOUT.inp +++ /dev/null @@ -1,39 +0,0 @@ -grid = fci.grid.nc - -nout = 1 -timestep = 0.01 - -MZ = 64 - -NXPE = 1 - -[mesh] -paralleltransform = fci - -symmetricglobalx = true - -[fci] -y_periodic = false -z_periodic = false - -[interpolation] -type=lagrange4pt - -[solver] -ATOL = 1e-12 -RTOL = 1e-8 -mms = true - -[f] -solution = sin(y - 2*z)*cos(t) + sin(y - z) - -source = -sin(t)*sin(y - 2*z) - (6.28318530717959*(0.01*x + 0.045)*(sin(y - z) + 2*cos(t)*cos(y - 2*z)) - 0.628318530717959*sin(y - z) - 0.628318530717959*cos(t)*cos(y - 2*z))/sqrt((0.01*x + 0.045)^2 + 1.0) - -bndry_par_all = parallel_dirichlet(f:solution) - -[g] -solution = -sin(y - 2*z)*cos(t) + cos(y - z) - -source = sin(t)*sin(y - 2*z) - (6.28318530717959*(0.01*x + 0.045)*(-2*cos(t)*cos(y - 2*z) - cos(y - z)) + 0.628318530717959*cos(t)*cos(y - 2*z) + 0.628318530717959*cos(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) - -bndry_par_all = parallel_dirichlet(g:solution) diff --git a/tests/integrated/test-fci-slab/plot_funcs.py b/tests/integrated/test-fci-slab/plot_funcs.py deleted file mode 100644 index d233858f8f..0000000000 --- a/tests/integrated/test-fci-slab/plot_funcs.py +++ /dev/null @@ -1,28 +0,0 @@ -from builtins import str -from builtins import range -# Plot interpolating functions -# Input generated by simple_test.py - -from numpy import linspace -import matplotlib.pyplot as plt -from boutdata.collect import collect - -f = collect("f", path="data") -yup = collect("yup", path="data") - -ny = 20 -nz = 8 - -# Note: yup[y=0] is never set -y = linspace(-1, 1, ny-1) - -plt.plot(f[0,4,4,:], 'o', label="f") - -for z in range(nz): - plt.plot(y+z, yup[4,1:,z], label="z = "+str(z)) - -plt.legend(loc='upper center') - -plt.savefig("plot_funcs.pdf") - -plt.show() diff --git a/tests/integrated/test-fci-slab/plot_interp.py b/tests/integrated/test-fci-slab/plot_interp.py deleted file mode 100644 index 2278684fba..0000000000 --- a/tests/integrated/test-fci-slab/plot_interp.py +++ /dev/null @@ -1,17 +0,0 @@ - -import matplotlib.pyplot as plt - -from boutdata.collect import collect - -f = collect("f", path="data") -yup = collect("yup", path="data") -ydown = collect("ydown", path="data") - -plt.plot(f[0,4,4,:], label="f") -plt.plot(yup[4,4,:], label="f.yup") -plt.plot(ydown[4,4,:], label="f.ydown") - -plt.legend() - -plt.savefig("plot_interp.pdf") -plt.show() diff --git a/tests/integrated/test-fci-slab/runtest b/tests/integrated/test-fci-slab/runtest deleted file mode 100755 index 196b2471f0..0000000000 --- a/tests/integrated/test-fci-slab/runtest +++ /dev/null @@ -1,165 +0,0 @@ -#!/usr/bin/env python3 -# -# Python script to run and analyse MMS test -# -from __future__ import division -from __future__ import print_function -from builtins import zip -from builtins import str - -from boututils.run_wrapper import shell, shell_safe, launch_safe, getmpirun -from boututils.datafile import DataFile -from boutdata.collect import collect - -from numpy import sqrt, max, abs, mean, array, log, pi, polyfit, linspace, arange - -import pickle - -from sys import stdout - -import zoidberg as zb - -showPlot = False #Do we want to show the plot as well as save it to file. - -nx = 5 # Not changed for these tests - -# Resolution in y and z -nlist = [64,128] #[8,16,32,64,128,256] - -nproc = 2 - -directory = "mms" - -varlist = ["f", "g"] -markers = ['bo', 'r^'] -labels = [r'$f$', r'$g$'] - -MPIRUN = getmpirun() - -success=True - -print("Making fci-slab test") -shell_safe("make > make.log") - -error_2 = {} -error_inf = {} -for var in varlist: - error_2[var] = [] # The L2 error (RMS) - error_inf[var] = [] # The maximum error - -yperiodic=False # Run with periodic Y? - -for n in nlist: - - # Define the magnetic field using new poloidal gridding method - # Note that the Bz and Bzprime parameters here must be the same as in mms.py - field = zb.field.Slab(Bz=0.05, Bzprime=0.1) - # Create rectangular poloidal grids - poloidal_grid = zb.poloidal_grid.RectangularPoloidalGrid(nx,n,1.,1.) - # Set the ylength and y locations - ylength = 10. - - if yperiodic: - ycoords = linspace(0.0, ylength, n, endpoint=False) - else: - # Doesn't include the end points - ycoords = (arange(n) + 0.5)*ylength/float(n) - - # Create the grid - grid = zb.grid.Grid(poloidal_grid, ycoords, ylength, yperiodic=yperiodic) - # Make and write maps - maps = zb.make_maps(grid, field) - zb.write_maps(grid, field, maps, new_names=False, metric2d=True) - - args = " -d "+directory+" MZ="+str(n)+ " fci:y_periodic="+str(yperiodic) - - # Command to run - cmd = "./fci_slab "+args - - print("Running command: "+cmd) - - # Launch using MPI - s, out = launch_safe(cmd, runcmd=MPIRUN, nproc=nproc, mthread=1, pipe=True) - - # Save output to log file - with open("run.log."+str(n), "w") as f: - f.write(out) - - if s: - print("Run failed!\nOutput was:\n") - print(out) - exit(s) - - for var in varlist: - # Collect data - E = collect("E_"+var, tind=[1,1], info=False, path=directory) - E = E[:,2:-2, :,:] - - # Average error over domain - l2 = sqrt(mean(E**2)) - linf = max(abs( E )) - - error_2[var].append( l2 ) - error_inf[var].append( linf ) - - print("%s : l-2 %f l-inf %f" % (var, l2, linf)) - -dx = 1. / array(nlist) - -# Save data -with open("fci_mms.pkl", "wb") as output: - pickle.dump(nlist, output) - pickle.dump(error_2, output) - pickle.dump(error_inf, output) - -# Calculate convergence order -for var,mark,label in zip(varlist, markers, labels): - fit = polyfit(log(dx), log(error_2[var]), 1) - order = fit[0] - stdout.write("%s Convergence order = %f (fit)" % (var, order)) - - order = log(error_2[var][-2]/error_2[var][-1])/log(dx[-2]/dx[-1]) - stdout.write(", %f (small spacing)" % (order,)) - - if order > 1.5: # Should be second order accurate - print("............ PASS") - else: - print("............ FAIL") - success = False - -if False: - try: - # Plot using matplotlib if available - import matplotlib.pyplot as plt - - plt.figure() - - for var,mark,label in zip(varlist, markers, labels): - plt.plot(dx, error_2[var], '-'+mark, label=label) - plt.plot(dx, error_inf[var], '--'+mark) - - plt.legend(loc="upper left") - plt.grid() - - plt.yscale('log') - plt.xscale('log') - - plt.xlabel(r'Mesh spacing $\delta x$') - plt.ylabel("Error norm") - - plt.savefig("fci-norm.pdf") - - print("Plot saved to fci-norm.pdf") - - if showPlot: - plt.show() - plt.close() - except ImportError: - print("No matplotlib") -else: - print("Plotting disabled") - -if success: - exit(0) -else: - exit(1) diff --git a/tests/integrated/test-fci-slab/simple_test.py b/tests/integrated/test-fci-slab/simple_test.py deleted file mode 100644 index 95f245c8c2..0000000000 --- a/tests/integrated/test-fci-slab/simple_test.py +++ /dev/null @@ -1,32 +0,0 @@ -from builtins import range -from numpy import zeros, linspace, concatenate -import boututils.datafile as bdata -from boutdata.input import transform3D - -# Parameters -nx = 10 -ny = 20 -nz = 8 - -shape = [nx, ny, nz] - -xt_prime = zeros(shape) -zt_prime = zeros(shape) - -for x in range(nx): - # No interpolation in x - xt_prime[x,:,:] = x - - # Each y slice scans between neighbouring z points - for z in range(nz): - zt_prime[x,:,z] = z + concatenate([linspace(-1, 1, ny-1), [0]]) - - -with bdata.DataFile('simple_test.nc', write=True, create=True) as f: - f.write('nx',nx) - f.write('ny',ny) - - for direction_name in ['forward', 'backward']: - f.write(direction_name + '_xt_prime', transform3D(xt_prime)) - f.write(direction_name + '_zt_prime', transform3D(zt_prime)) - diff --git a/tests/integrated/test-smooth/test_smooth.cxx b/tests/integrated/test-smooth/test_smooth.cxx index f1246d0fed..c4670656c1 100644 --- a/tests/integrated/test-smooth/test_smooth.cxx +++ b/tests/integrated/test-smooth/test_smooth.cxx @@ -17,7 +17,7 @@ int main(int argc, char **argv) { Field2D input2d = f.create2D("1 + sin(2*y)"); Field3D input3d = f.create3D("gauss(x-0.5,0.2)*gauss(y-pi)*sin(3*y - z)"); - input3d.mergeYupYdown(); + mesh->getParallelTransform().calcYUpDown(input3d); SAVE_ONCE2(input2d, input3d); @@ -31,11 +31,6 @@ int main(int argc, char **argv) { // Output data dump.write(); - dump.close(); - - output << "\nFinished running test. Triggering error to quit\n\n"; - - MPI_Barrier(BoutComm::get()); // Wait for all processors to write data BoutFinalise(); return 0; diff --git a/tests/unit/field/test_field3d.cxx b/tests/unit/field/test_field3d.cxx index 7f62617c36..f7260739c8 100644 --- a/tests/unit/field/test_field3d.cxx +++ b/tests/unit/field/test_field3d.cxx @@ -223,20 +223,15 @@ TEST_F(Field3DTest, MergeYupYDown) { field.mergeYupYdown(); - EXPECT_TRUE(field.hasYupYdown()); + EXPECT_FALSE(field.hasYupYdown()); - auto& yup = field.yup(); - EXPECT_EQ(&field, &yup); - auto& ydown = field.ydown(); - EXPECT_EQ(&field, &ydown); +#if CHECK > 2 + EXPECT_THROW(field.yup(), BoutException); + EXPECT_THROW(field.ydown(), BoutException); +#endif // Should be able to merge again without any problems - field.mergeYupYdown(); - - auto& yup2 = field.yup(); - EXPECT_EQ(&field, &yup2); - auto& ydown2 = field.ydown(); - EXPECT_EQ(&field, &ydown2); + EXPECT_NO_THROW(field.mergeYupYdown()); } TEST_F(Field3DTest, SplitThenMergeYupYDown) { @@ -252,10 +247,36 @@ TEST_F(Field3DTest, SplitThenMergeYupYDown) { field.mergeYupYdown(); - auto& yup2 = field.yup(); - EXPECT_EQ(&field, &yup2); - auto& ydown2 = field.ydown(); - EXPECT_EQ(&field, &ydown2); +#if CHECK > 2 + EXPECT_THROW(field.yup(), BoutException); + EXPECT_THROW(field.ydown(), BoutException); +#endif +} + +TEST_F(Field3DTest, MultipleYupYdown) { + FakeMesh newmesh{3, 5, 7}; + newmesh.ystart = 2; + + Field3D field{&newmesh}; + + field.splitYupYdown(); + + EXPECT_TRUE(field.hasYupYdown()); + + auto &yup = field.yup(); + EXPECT_NE(&field, &yup); + auto &ydown = field.ydown(); + EXPECT_NE(&field, &ydown); + auto &yup1 = field.yup(1); + EXPECT_NE(&field, &yup1); + EXPECT_NE(&yup, &yup1); + auto &ydown1 = field.ydown(1); + EXPECT_NE(&field, &ydown1); + EXPECT_NE(&ydown, &ydown1); + +#if CHECK > 1 + EXPECT_THROW(field.yup(2), BoutException); +#endif } TEST_F(Field3DTest, Ynext) { @@ -270,7 +291,9 @@ TEST_F(Field3DTest, Ynext) { EXPECT_NE(&field, &ydown); EXPECT_NE(&yup, &ydown); +#if CHECK > 0 EXPECT_THROW(field.ynext(99), BoutException); +#endif } TEST_F(Field3DTest, ConstYnext) { @@ -286,7 +309,9 @@ TEST_F(Field3DTest, ConstYnext) { EXPECT_NE(&field2, &ydown); EXPECT_NE(&yup, &ydown); +#if CHECK > 0 EXPECT_THROW(field2.ynext(99), BoutException); +#endif } TEST_F(Field3DTest, GetGlobalMesh) { diff --git a/tests/unit/invert/test_fft.cxx b/tests/unit/invert/test_fft.cxx index 37ed3cb9b8..735fcedb93 100644 --- a/tests/unit/invert/test_fft.cxx +++ b/tests/unit/invert/test_fft.cxx @@ -39,9 +39,6 @@ class FFTTest : public ::testing::TestWithParam { Array real_signal; Array fft_signal; - - // FFTs have a slightly looser tolerance than other functions - static constexpr BoutReal fft_tolerance{1.e-12}; }; // Test the FFT functions with both even- and odd-length real signals @@ -57,8 +54,8 @@ TEST_P(FFTTest, rfft) { EXPECT_EQ(output.size(), nmodes); for (int i = 0; i < nmodes; ++i) { - EXPECT_NEAR(real(output[i]), real(fft_signal[i]), fft_tolerance); - EXPECT_NEAR(imag(output[i]), imag(fft_signal[i]), fft_tolerance); + EXPECT_NEAR(real(output[i]), real(fft_signal[i]), FFTTolerance); + EXPECT_NEAR(imag(output[i]), imag(fft_signal[i]), FFTTolerance); } } @@ -72,7 +69,7 @@ TEST_P(FFTTest, irfft) { EXPECT_EQ(output.size(), size); for (int i = 0; i < size; ++i) { - EXPECT_NEAR(output[i], real_signal[i], fft_tolerance); + EXPECT_NEAR(output[i], real_signal[i], FFTTolerance); } } @@ -84,8 +81,8 @@ TEST_P(FFTTest, rfftWithArray) { EXPECT_EQ(output.size(), nmodes); for (int i = 0; i < nmodes; ++i) { - EXPECT_NEAR(real(output[i]), real(fft_signal[i]), fft_tolerance); - EXPECT_NEAR(imag(output[i]), imag(fft_signal[i]), fft_tolerance); + EXPECT_NEAR(real(output[i]), real(fft_signal[i]), FFTTolerance); + EXPECT_NEAR(imag(output[i]), imag(fft_signal[i]), FFTTolerance); } } @@ -97,7 +94,7 @@ TEST_P(FFTTest, irfftWithArray) { EXPECT_EQ(output.size(), size); for (int i = 0; i < size; ++i) { - EXPECT_NEAR(output[i], real_signal[i], fft_tolerance); + EXPECT_NEAR(output[i], real_signal[i], FFTTolerance); } } @@ -109,6 +106,6 @@ TEST_P(FFTTest, RoundTrip) { EXPECT_EQ(output.size(), real_signal.size()); for (int i = 0; i < size; ++i) { - EXPECT_NEAR(output[i], real_signal[i], fft_tolerance); + EXPECT_NEAR(output[i], real_signal[i], FFTTolerance); } } diff --git a/tests/unit/mesh/parallel/test_shiftedmetric.cxx b/tests/unit/mesh/parallel/test_shiftedmetric.cxx index aed252d2f4..37d13f6179 100644 --- a/tests/unit/mesh/parallel/test_shiftedmetric.cxx +++ b/tests/unit/mesh/parallel/test_shiftedmetric.cxx @@ -1,7 +1,8 @@ #include "gtest/gtest.h" -#include "bout/paralleltransform.hxx" +#include "fft.hxx" #include "test_extras.hxx" +#include "bout/paralleltransform.hxx" // The unit tests use the global mesh using namespace bout::globals; @@ -25,9 +26,42 @@ class ShiftedMetricTest : public ::testing::Test { mesh->createDefaultRegions(); output_info.enable(); + // Make sure fft functions are quiet by setting fft_measure to false + bout::fft::fft_init(false); + zShift = Field2D{mesh}; - fillField(zShift, {{1., 2., 3., 4., 5.}, {1., 2., 3., 4., 5.}, {1., 2., 3., 4., 5.}}); + fillField(zShift, {{1., 2., 3., 4., 5., 6., 7.}, + {2., 4., 6., 8., 10., 12., 14.}, + {3., 6., 9., 12., 15., 18., 21.}}); + + Field3D input_temp{mesh}; + + fillField(input_temp, {{{1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}}, + + {{1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}}, + + {{1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}}}); + + input = std::move(input_temp); dynamic_cast(mesh)->setCoordinates(std::make_shared( mesh, Field2D{1.0}, Field2D{1.0}, BoutReal{1.0}, Field2D{1.0}, Field2D{0.0}, @@ -42,182 +76,221 @@ class ShiftedMetricTest : public ::testing::Test { } static constexpr int nx = 3; - static constexpr int ny = 5; - static constexpr int nz = 7; + static constexpr int ny = 7; + static constexpr int nz = 5; Field2D zShift; + Field3D input; }; TEST_F(ShiftedMetricTest, ToFieldAligned) { ShiftedMetric shifted{*mesh, zShift}; - Field3D input{mesh}; - - fillField(input, {{{1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}}, - - {{1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}}, - - {{1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}}}); - Field3D expected{mesh}; - fillField(expected, {{{2., 3., 4., 5., 6., 7., 1.}, - {3., 4., 5., 6., 7., 1., 2.}, - {4., 5., 6., 7., 1., 2., 3.}, - {5., 6., 7., 1., 2., 3., 4.}, - {6., 7., 1., 2., 3., 4., 5.}}, - - {{2., 3., 4., 5., 6., 7., 1.}, - {3., 4., 5., 6., 7., 1., 2.}, - {4., 5., 6., 7., 1., 2., 3.}, - {5., 6., 7., 1., 2., 3., 4.}, - {6., 7., 1., 2., 3., 4., 5.}}, - - {{2., 3., 4., 5., 6., 7., 1.}, - {3., 4., 5., 6., 7., 1., 2.}, - {4., 5., 6., 7., 1., 2., 3.}, - {5., 6., 7., 1., 2., 3., 4.}, - {6., 7., 1., 2., 3., 4., 5.}}}); - - EXPECT_TRUE(IsFieldEqual(shifted.toFieldAligned(input), expected)); + fillField(expected, {{{2., 3., 4., 5., 1.}, + {3., 4., 5., 1., 2.}, + {4., 5., 1., 2., 3.}, + {5., 1., 2., 3., 4.}, + {1., 2., 3., 4., 5.}, + {2., 3., 4., 5., 1.}, + {3., 4., 5., 1., 2.}}, + + {{3., 4., 5., 1., 2.}, + {5., 1., 2., 3., 4.}, + {2., 3., 4., 5., 1.}, + {4., 5., 1., 2., 3.}, + {1., 2., 3., 4., 5.}, + {3., 4., 5., 1., 2.}, + {5., 1., 2., 3., 4.}}, + + {{4., 5., 1., 2., 3.}, + {2., 3., 4., 5., 1.}, + {5., 1., 2., 3., 4.}, + {3., 4., 5., 1., 2.}, + {1., 2., 3., 4., 5.}, + {4., 5., 1., 2., 3.}, + {2., 3., 4., 5., 1.}}}); + + EXPECT_TRUE(IsFieldEqual(shifted.toFieldAligned(input), expected, "RGN_ALL", + FFTTolerance)); } TEST_F(ShiftedMetricTest, FromFieldAligned) { ShiftedMetric shifted{*mesh, zShift}; - Field3D input{mesh}; - - fillField(input, {{{1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}}, - - {{1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}}, - - {{1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}}}); - Field3D expected{mesh}; - fillField(expected, {{{7., 1., 2., 3., 4., 5., 6.}, - {6., 7., 1., 2., 3., 4., 5.}, - {5., 6., 7., 1., 2., 3., 4.}, - {4., 5., 6., 7., 1., 2., 3.}, - {3., 4., 5., 6., 7., 1., 2.}}, - - {{7., 1., 2., 3., 4., 5., 6.}, - {6., 7., 1., 2., 3., 4., 5.}, - {5., 6., 7., 1., 2., 3., 4.}, - {4., 5., 6., 7., 1., 2., 3.}, - {3., 4., 5., 6., 7., 1., 2.}}, - - {{7., 1., 2., 3., 4., 5., 6.}, - {6., 7., 1., 2., 3., 4., 5.}, - {5., 6., 7., 1., 2., 3., 4.}, - {4., 5., 6., 7., 1., 2., 3.}, - {3., 4., 5., 6., 7., 1., 2.}}}); + fillField(expected, {{{5., 1., 2., 3., 4.}, + {4., 5., 1., 2., 3.}, + {3., 4., 5., 1., 2.}, + {2., 3., 4., 5., 1.}, + {1., 2., 3., 4., 5.}, + {5., 1., 2., 3., 4.}, + {4., 5., 1., 2., 3.}}, + + {{4., 5., 1., 2., 3.}, + {2., 3., 4., 5., 1.}, + {5., 1., 2., 3., 4.}, + {3., 4., 5., 1., 2.}, + {1., 2., 3., 4., 5.}, + {4., 5., 1., 2., 3.}, + {2., 3., 4., 5., 1.}}, + + {{3., 4., 5., 1., 2.}, + {5., 1., 2., 3., 4.}, + {2., 3., 4., 5., 1.}, + {4., 5., 1., 2., 3.}, + {1., 2., 3., 4., 5.}, + {3., 4., 5., 1., 2.}, + {5., 1., 2., 3., 4.}}}); // Loosen tolerance a bit due to FFTs - EXPECT_TRUE(IsFieldEqual(shifted.fromFieldAligned(input), expected, "RGN_ALL", 1.e-12)); + EXPECT_TRUE(IsFieldEqual(shifted.fromFieldAligned(input), expected, "RGN_ALL", + FFTTolerance)); } TEST_F(ShiftedMetricTest, CalcYUpDown) { + // Use two y-guards to test multiple parallel slices + mesh->ystart = 2; + mesh->yend = mesh->LocalNy - 3; + + // We don't shift in the guard cells, and the parallel slices are + // stored offset in y, therefore we need to make new regions that we + // can compare the expected and actual outputs over output_info.disable(); - auto region_yup = mesh->getRegion("RGN_NOY"); - region_yup.periodicShift(ShiftedMetricTest::nz, - ShiftedMetricTest::ny * ShiftedMetricTest::nz); - mesh->addRegion("RGN_YUP", region_yup); - - auto region_ydown = mesh->getRegion("RGN_NOY"); - region_ydown.periodicShift(-ShiftedMetricTest::nz, - ShiftedMetricTest::ny * ShiftedMetricTest::nz); - mesh->addRegion("RGN_YDOWN", region_ydown); + mesh->addRegion3D("RGN_YUP", + Region(0, mesh->LocalNx - 1, mesh->ystart + 1, mesh->yend + 1, + 0, mesh->LocalNz - 1, mesh->LocalNy, mesh->LocalNz)); + mesh->addRegion3D("RGN_YUP2", + Region(0, mesh->LocalNx - 1, mesh->ystart + 2, mesh->yend + 2, + 0, mesh->LocalNz - 1, mesh->LocalNy, mesh->LocalNz)); + + mesh->addRegion3D("RGN_YDOWN", + Region(0, mesh->LocalNx - 1, mesh->ystart - 1, mesh->yend - 1, + 0, mesh->LocalNz - 1, mesh->LocalNy, mesh->LocalNz)); + mesh->addRegion3D("RGN_YDOWN2", + Region(0, mesh->LocalNx - 1, mesh->ystart - 2, mesh->yend - 2, + 0, mesh->LocalNz - 1, mesh->LocalNy, mesh->LocalNz)); output_info.enable(); + // Actual interesting bit here! ShiftedMetric shifted{*mesh, zShift}; - - Field3D input{mesh}; - - fillField(input, {{{1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}}, - - {{1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}}, - - {{1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}}}); - shifted.calcYUpDown(input); - Field3D expected_up{mesh}; - - fillField(expected_up, {{{0., 0., 0., 0., 0., 0., 0.}, - {2., 3., 4., 5., 6., 7., 1.}, - {2., 3., 4., 5., 6., 7., 1.}, - {2., 3., 4., 5., 6., 7., 1.}, - {2., 3., 4., 5., 6., 7., 1.}}, - - {{0., 0., 0., 0., 0., 0., 0.}, - {2., 3., 4., 5., 6., 7., 1.}, - {2., 3., 4., 5., 6., 7., 1.}, - {2., 3., 4., 5., 6., 7., 1.}, - {2., 3., 4., 5., 6., 7., 1.}}, - - {{0., 0., 0., 0., 0., 0., 0.}, - {2., 3., 4., 5., 6., 7., 1.}, - {2., 3., 4., 5., 6., 7., 1.}, - {2., 3., 4., 5., 6., 7., 1.}, - {2., 3., 4., 5., 6., 7., 1.}}}); - - Field3D expected_down{mesh}; - - fillField(expected_down, {{{7., 1., 2., 3., 4., 5., 6.}, - {7., 1., 2., 3., 4., 5., 6.}, - {7., 1., 2., 3., 4., 5., 6.}, - {7., 1., 2., 3., 4., 5., 6.}, - {0., 0., 0., 0., 0., 0., 0.}}, - - {{7., 1., 2., 3., 4., 5., 6.}, - {7., 1., 2., 3., 4., 5., 6.}, - {7., 1., 2., 3., 4., 5., 6.}, - {7., 1., 2., 3., 4., 5., 6.}, - {0., 0., 0., 0., 0., 0., 0.}}, - - {{7., 1., 2., 3., 4., 5., 6.}, - {7., 1., 2., 3., 4., 5., 6.}, - {7., 1., 2., 3., 4., 5., 6.}, - {7., 1., 2., 3., 4., 5., 6.}, - {0., 0., 0., 0., 0., 0., 0.}}}); - - EXPECT_TRUE(IsFieldEqual(input.yup(), expected_up, "RGN_YUP")); - EXPECT_TRUE(IsFieldEqual(input.ydown(), expected_down, "RGN_YDOWN")); + // Expected output values + + Field3D expected_up_1{mesh}; + + // Note: here zeroes are for values we don't expect to read + fillField(expected_up_1, {{{0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {2., 3., 4., 5., 1.}, + {2., 3., 4., 5., 1.}, + {2., 3., 4., 5., 1.}, + {0., 0., 0., 0., 0.}}, + + {{0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {3., 4., 5., 1., 2.}, + {3., 4., 5., 1., 2.}, + {3., 4., 5., 1., 2.}, + {0., 0., 0., 0., 0.}}, + + {{0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {4., 5., 1., 2., 3.}, + {4., 5., 1., 2., 3.}, + {4., 5., 1., 2., 3.}, + {0., 0., 0., 0., 0.}}}); + + Field3D expected_up_2{mesh}; + + fillField(expected_up_2, {{{0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {3., 4., 5., 1., 2.}, + {3., 4., 5., 1., 2.}, + {3., 4., 5., 1., 2.}}, + + {{0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {5., 1., 2., 3., 4.}, + {5., 1., 2., 3., 4.}, + {5., 1., 2., 3., 4.}}, + + {{0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {2., 3., 4., 5., 1.}, + {2., 3., 4., 5., 1.}, + {2., 3., 4., 5., 1.}}}); + + Field3D expected_down_1{mesh}; + + fillField(expected_down_1, {{{0., 0., 0., 0., 0.}, + {5., 1., 2., 3., 4.}, + {5., 1., 2., 3., 4.}, + {5., 1., 2., 3., 4.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}}, + + {{0., 0., 0., 0., 0.}, + {4., 5., 1., 2., 3.}, + {4., 5., 1., 2., 3.}, + {4., 5., 1., 2., 3.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}}, + + {{0., 0., 0., 0., 0.}, + {3., 4., 5., 1., 2.}, + {3., 4., 5., 1., 2.}, + {3., 4., 5., 1., 2.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}}}); + + Field3D expected_down2{mesh}; + + fillField(expected_down2, {{{4., 5., 1., 2., 3.}, + {4., 5., 1., 2., 3.}, + {4., 5., 1., 2., 3.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}}, + + {{2., 3., 4., 5., 1.}, + {2., 3., 4., 5., 1.}, + {2., 3., 4., 5., 1.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}}, + + {{5., 1., 2., 3., 4.}, + {5., 1., 2., 3., 4.}, + {5., 1., 2., 3., 4.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}}}); + + EXPECT_TRUE( + IsFieldEqual(input.ynext(1), expected_up_1, "RGN_YUP", FFTTolerance)); + EXPECT_TRUE( + IsFieldEqual(input.ynext(2), expected_up_2, "RGN_YUP2", FFTTolerance)); + EXPECT_TRUE( + IsFieldEqual(input.ynext(-1), expected_down_1, "RGN_YDOWN", FFTTolerance)); + EXPECT_TRUE( + IsFieldEqual(input.ynext(-2), expected_down2, "RGN_YDOWN2", FFTTolerance)); } diff --git a/tests/unit/mesh/test_paralleltransform.cxx b/tests/unit/mesh/test_paralleltransform.cxx new file mode 100644 index 0000000000..3996db516f --- /dev/null +++ b/tests/unit/mesh/test_paralleltransform.cxx @@ -0,0 +1,41 @@ +#include "gtest/gtest.h" + +#include "test_extras.hxx" +#include "bout/paralleltransform.hxx" + +namespace bout { +namespace globals { +extern Mesh* mesh; +} +} // namespace bout + +using ParallelTransformTest = FakeMeshFixture; + +TEST_F(ParallelTransformTest, IdentityCalcYUpDown) { + + ParallelTransformIdentity transform{}; + + Field3D field{1.0}; + + transform.calcYUpDown(field); + + EXPECT_TRUE(IsFieldEqual(field.yup(), 1.0)); + EXPECT_TRUE(IsFieldEqual(field.ydown(), 1.0)); +} + +TEST_F(ParallelTransformTest, IdentityCalcYUpDownTwoSlices) { + + ParallelTransformIdentity transform{}; + + bout::globals::mesh->ystart = 2; + + Field3D field{1.0}; + + transform.calcYUpDown(field); + + EXPECT_TRUE(IsFieldEqual(field.yup(0), 1.0)); + EXPECT_TRUE(IsFieldEqual(field.yup(1), 1.0)); + + EXPECT_TRUE(IsFieldEqual(field.ydown(0), 1.0)); + EXPECT_TRUE(IsFieldEqual(field.ydown(1), 1.0)); +} diff --git a/tests/unit/test_extras.hxx b/tests/unit/test_extras.hxx index 4c483fc037..734edd6e10 100644 --- a/tests/unit/test_extras.hxx +++ b/tests/unit/test_extras.hxx @@ -13,7 +13,9 @@ #include "field3d.hxx" #include "unused.hxx" -const BoutReal BoutRealTolerance = 1e-15; +static constexpr BoutReal BoutRealTolerance{1e-15}; +// FFTs have a slightly looser tolerance than other functions +static constexpr BoutReal FFTTolerance{1.e-12}; /// Does \p str contain \p substring? ::testing::AssertionResult IsSubString(const std::string &str, diff --git a/tools/pylib/zoidberg/grid.py b/tools/pylib/zoidberg/grid.py index 99e247d43f..8ca964517b 100644 --- a/tools/pylib/zoidberg/grid.py +++ b/tools/pylib/zoidberg/grid.py @@ -184,8 +184,8 @@ def metric(self): # Note: These y metrics are for Cartesian coordinates # If in cylindrical coordinates then these should be different - g_yy = 1.0 # Rmaj**2 - gyy = 1.0 # 1/Rmaj**2 + g_yy = np.ones(self.shape) + gyy = np.ones(self.shape) return {"dx":dx, "dy":dy3d, "dz": dz, "gyy": gyy, "g_yy":g_yy, diff --git a/tools/pylib/zoidberg/test_zoidberg.py b/tools/pylib/zoidberg/test_zoidberg.py index d82c30f0e2..45c1626381 100644 --- a/tools/pylib/zoidberg/test_zoidberg.py +++ b/tools/pylib/zoidberg/test_zoidberg.py @@ -1,69 +1,68 @@ - +from itertools import chain, product import numpy as np from . import zoidberg, grid, field + def test_make_maps_slab(): nx = 5 - ny = 6 + ny = 8 nz = 7 # Create a straight magnetic field in a slab straight_field = field.Slab(By=1.0, Bz=0.0, Bzprime=0.0) - + # Create a rectangular grid in (x,y,z) - rectangle = grid.rectangular_grid(nx,ny,nz) + rectangle = grid.rectangular_grid(nx, ny, nz) + + # Two parallel slices in each direction + nslice = 2 # Calculate forwards and backwards maps - maps = zoidberg.make_maps(rectangle, straight_field) - - # Check that maps has the required forward and backward index variables - for var in ['forward_xt_prime', 'forward_zt_prime', 'backward_xt_prime', 'backward_zt_prime']: - assert var in maps + maps = zoidberg.make_maps(rectangle, straight_field, nslice=nslice) - # Each map should have the same shape as the grid - assert maps['forward_xt_prime'].shape == (nx,ny,nz) - assert maps['backward_xt_prime'].shape == (nx,ny,nz) - assert maps['forward_zt_prime'].shape == (nx,ny,nz) - assert maps['backward_zt_prime'].shape == (nx,ny,nz) - - # Since this is a straight magnetic field in a simple rectangle, + # Since this is a straight magnetic field in a simple rectangle, # all the maps should be the same, and should be the identity - identity_map_x, identity_map_z = np.meshgrid(np.arange(nx), np.arange(nz), indexing='ij') - for y in range(ny-1): - assert np.allclose(maps['forward_xt_prime'][:,y,:], identity_map_x) - assert np.allclose(maps['forward_zt_prime'][:,y,:], identity_map_z) + # Check that maps has the required forward and backward index variables + offsets = chain(range(1, nslice + 1), range(-1, -(nslice + 1), -1)) + field_line_maps = ["xt_prime", "zt_prime"] + + for field_line_map, offset in product(field_line_maps, offsets): + var = zoidberg.parallel_slice_field_name(field_line_map, offset) + print("Current field: ", var) + assert var in maps + + # Each map should have the same shape as the grid + assert maps[var].shape == (nx, ny, nz) - for y in range(1,ny): - assert np.allclose(maps['backward_xt_prime'][:,y,:], identity_map_x) - assert np.allclose(maps['backward_zt_prime'][:,y,:], identity_map_z) + # The first/last abs(offset) points are not valid, so ignore those + interior_range = range(ny-abs(offset)) if offset > 0 else range(abs(offset), ny) + # Those invalid points should be set to -1 + end_slice = slice(-1, -(offset + 1), -1) if offset > 0 else slice(0, -offset) + identity_map = identity_map_x if "x" in var else identity_map_z - # The last forward map should hit a boundary - assert np.allclose(maps['forward_xt_prime'][:,-1,:], -1.0) - assert np.allclose(maps['forward_zt_prime'][:,-1,:], -1.0) + for y in interior_range: + assert np.allclose(maps[var][:, y, :], identity_map) - # First backward map hits boundary - assert np.allclose(maps['backward_xt_prime'][:,0,:], -1.0) - assert np.allclose(maps['backward_zt_prime'][:,0,:], -1.0) + # The end slice should hit a boundary + assert np.allclose(maps[var][:, end_slice, :], -1.0) def test_make_maps_straight_stellarator(): nx = 5 ny = 6 nz = 7 - + # Create magnetic field magnetic_field = field.StraightStellarator(radius = np.sqrt(2.0)) - + # Create a rectangular grid in (x,y,z) rectangle = grid.rectangular_grid(nx,ny,nz, Lx = 1.0, Lz = 1.0, Ly = 10.0, yperiodic = True) - + # Here both the field and and grid are centred at (x,z) = (0,0) # and the rectangular grid here fits entirely within the coils maps = zoidberg.make_maps(rectangle, magnetic_field) - - diff --git a/tools/pylib/zoidberg/zoidberg.py b/tools/pylib/zoidberg/zoidberg.py index a7333d6762..7d11882fb2 100644 --- a/tools/pylib/zoidberg/zoidberg.py +++ b/tools/pylib/zoidberg/zoidberg.py @@ -2,6 +2,11 @@ import numpy as np from boututils import datafile as bdata +from collections import namedtuple +from itertools import chain + +from . import fieldtracer +from .progress import update_progress # PyEVTK might be called pyevtk or evtk, depending on where it was # installed from @@ -14,13 +19,24 @@ except ImportError: have_evtk = False -# from . import grid -# from . import field -from . import fieldtracer -from .progress import update_progress + +def parallel_slice_field_name(field, offset): + """Form a unique, backwards-compatible name for field at a given offset + + Parameters + ---------- + field : str + Name of the field to convert + offset : int + Parallel slice offset + + """ + prefix = 'forward' if offset > 0 else 'backward' + suffix = "_{}".format(abs(offset)) if abs(offset) > 1 else "" + return "{}_{}{}".format(prefix, field, suffix) -def make_maps(grid, magnetic_field, quiet=False, **kwargs): +def make_maps(grid, magnetic_field, nslice=1, quiet=False, **kwargs): """Make the forward and backward FCI maps Parameters @@ -29,6 +45,8 @@ def make_maps(grid, magnetic_field, quiet=False, **kwargs): Grid generated by Zoidberg magnetic_field : :py:obj:`zoidberg.field.MagneticField` Zoidberg magnetic field object + nslice : int + Number of parallel slices in each direction quiet : bool Don't display progress bar kwargs @@ -51,115 +69,91 @@ def make_maps(grid, magnetic_field, quiet=False, **kwargs): shape = (nx, ny, nz) # Coordinates of each grid point - R = np.zeros( shape ) - Z = np.zeros( shape ) - - # Arrays to store X index at end of field-line - # starting from (x,y,z) and going forward in toroidal angle (y) - forward_xt_prime = np.zeros( shape ) - forward_zt_prime = np.zeros( shape ) - - forward_R = np.zeros( shape ) - forward_Z = np.zeros( shape ) - - # Same but going backwards in toroidal angle - backward_xt_prime = np.zeros( shape ) - backward_zt_prime = np.zeros( shape ) - - backward_R = np.zeros( shape ) - backward_Z = np.zeros( shape ) - - field_tracer = fieldtracer.FieldTracer(magnetic_field) - - try: - rtol = kwargs["rtol"] - except KeyError: - rtol = None + R = np.zeros(shape) + Z = np.zeros(shape) - # TODO: if axisymmetric, don't loop, do one slice and copy for j in range(ny): - if (not quiet) and (ny > 1): - update_progress(float(j)/float(ny-1), **kwargs) - - # Get this poloidal grid - pol, ycoord = grid.getPoloidalGrid(j) - - # Store coordinates - R[:,j,:] = pol.R - Z[:,j,:] = pol.Z - - # Get the next (forward) poloidal grid - pol_forward, y_forward = grid.getPoloidalGrid(j+1) - - # We only want the end point, as [0,...] is the initial position - coord = field_tracer.follow_field_lines(pol.R, pol.Z, [ycoord, y_forward], rtol=rtol)[1,...] - - # Store the coordinates in real space - forward_R[:,j,:] = coord[:,:,0] - forward_Z[:,j,:] = coord[:,:,1] - - # Get the indices into the forward poloidal grid - if pol_forward is None: - # No forward grid, so hit a boundary - xind = -1 - zind = -1 - else: - # Find the indices for these new locations on the forward poloidal grid - xcoord = coord[:,:,0] - zcoord = coord[:,:,1] - xind, zind = pol_forward.findIndex(xcoord, zcoord) - - # Check boundary defined by the field - outside = magnetic_field.boundary.outside(xcoord, y_forward, zcoord) - xind[outside] = -1 - zind[outside] = -1 - - forward_xt_prime[:,j,:] = xind - forward_zt_prime[:,j,:] = zind - - # Go backwards one poloidal grid - pol_back, y_back = grid.getPoloidalGrid(j-1) - - # We only want the end point, as [0,...] is the initial position - coord = field_tracer.follow_field_lines(pol.R, pol.Z, [ycoord, y_back], rtol=rtol)[1,...] - - # Store the coordinates in real space - backward_R[:,j,:] = coord[:,:,0] - backward_Z[:,j,:] = coord[:,:,1] - - if pol_back is None: - # Hit boundary - xind = -1 - zind = -1 - else: - # Find the indices for these new locations on the backward poloidal grid - xcoord = coord[:,:,0] - zcoord = coord[:,:,1] - xind, zind = pol_back.findIndex(xcoord, zcoord) - - # Check boundary defined by the field - outside = magnetic_field.boundary.outside(xcoord, y_back, zcoord) - xind[outside] = -1 - zind[outside] = -1 - - backward_xt_prime[:,j,:] = xind - backward_zt_prime[:,j,:] = zind + pol, _ = grid.getPoloidalGrid(j) + R[:, j, :] = pol.R + Z[:, j, :] = pol.Z + + field_tracer = fieldtracer.FieldTracer(magnetic_field) + + rtol = kwargs.get("rtol", None) + # The field line maps and coordinates, etc. maps = { - 'R' : R, 'Z':Z, - 'forward_R':forward_R, 'forward_Z':forward_Z, - 'forward_xt_prime' : forward_xt_prime, - 'forward_zt_prime' : forward_zt_prime, - 'backward_R':backward_R, 'backward_Z':backward_Z, - 'backward_xt_prime' : backward_xt_prime, - 'backward_zt_prime' : backward_zt_prime + 'R': R, + 'Z': Z, } + # A helper data structure that groups the various field line maps along with the offset + ParallelSlice = namedtuple('ParallelSlice', ['offset', 'R', 'Z', 'xt_prime', 'zt_prime']) + # A list of the above data structures for each offset we want + parallel_slices = [] + + # Loop over offsets {1, ... nslice, -1, ... -nslice} + for offset in chain(range(1, nslice + 1), range(-1, -(nslice + 1), -1)): + # Unique names of the field line maps for this offset + field_names = [parallel_slice_field_name(field, offset) + for field in ['R', 'Z', 'xt_prime', 'zt_prime']] + + # Initialise the field arrays -- puts them straight into the result dict + for field in field_names: + maps[field] = np.zeros(shape) + + # Get the field arrays we just made and wrap them up in our helper tuple + fields = map(lambda x: maps[x], field_names) + parallel_slices.append(ParallelSlice(offset, *fields)) + + # Total size of the progress bar + total_work = float((len(parallel_slices) - 1) * (ny-1)) + + # TODO: if axisymmetric, don't loop, do one slice and copy + # TODO: restart tracing for adjacent offsets + for slice_index, parallel_slice in enumerate(parallel_slices): + for j in range(ny): + if (not quiet) and (ny > 1): + update_progress(float(slice_index * j) / total_work, **kwargs) + + # Get this poloidal grid + pol, ycoord = grid.getPoloidalGrid(j) + + # Get the next poloidal grid + pol_slice, y_slice = grid.getPoloidalGrid(j + parallel_slice.offset) + + # We only want the end point, as [0,...] is the initial position + coord = field_tracer.follow_field_lines(pol.R, pol.Z, [ycoord, y_slice], rtol=rtol)[1, ...] + + # Store the coordinates in real space + parallel_slice.R[:, j, :] = coord[:, :, 0] + parallel_slice.Z[:, j, :] = coord[:, :, 1] + + # Get the indices into the slice poloidal grid + if pol_slice is None: + # No slice grid, so hit a boundary + xind = -1 + zind = -1 + else: + # Find the indices for these new locations on the slice poloidal grid + xcoord = coord[:, :, 0] + zcoord = coord[:, :, 1] + xind, zind = pol_slice.findIndex(xcoord, zcoord) + + # Check boundary defined by the field + outside = magnetic_field.boundary.outside(xcoord, y_slice, zcoord) + xind[outside] = -1 + zind[outside] = -1 + + parallel_slice.xt_prime[:, j, :] = xind + parallel_slice.zt_prime[:, j, :] = zind + return maps def write_maps(grid, magnetic_field, maps, gridfile='fci.grid.nc', - new_names=False, metric2d=True, format="NETCDF3_64BIT"): + new_names=False, metric2d=True, format="NETCDF3_64BIT", + quiet=False): """Write FCI maps to BOUT++ grid file Parameters @@ -176,8 +170,10 @@ def write_maps(grid, magnetic_field, maps, gridfile='fci.grid.nc', Write "g_yy" rather than "g_22" metric2d : bool, optional Output only 2D metrics - format : str + format : str, optional Specifies file format to use, passed to boutdata.DataFile + quiet : bool, optional + Don't warn about 2D metrics Returns ------- @@ -193,24 +189,31 @@ def write_maps(grid, magnetic_field, maps, gridfile='fci.grid.nc', # Check if the magnetic field is in cylindrical coordinates # If so, we need to change the gyy and g_yy metrics - pol_grid,ypos = grid.getPoloidalGrid(0) + pol_grid, ypos = grid.getPoloidalGrid(0) Rmaj = magnetic_field.Rfunc(pol_grid.R, pol_grid.Z, ypos) if Rmaj is not None: # In cylindrical coordinates Rmaj = np.zeros(grid.shape) for yindex in range(grid.numberOfPoloidalGrids()): - pol_grid,ypos = grid.getPoloidalGrid(yindex) - Rmaj[:,yindex,:] = magnetic_field.Rfunc(pol_grid.R, pol_grid.Z, ypos) + pol_grid, ypos = grid.getPoloidalGrid(yindex) + Rmaj[:, yindex, :] = magnetic_field.Rfunc(pol_grid.R, pol_grid.Z, ypos) metric["gyy"] = 1./Rmaj**2 metric["g_yy"] = Rmaj**2 - + # Get magnetic field and pressure Bmag = np.zeros(grid.shape) pressure = np.zeros(grid.shape) for yindex in range(grid.numberOfPoloidalGrids()): - pol_grid,ypos = grid.getPoloidalGrid(yindex) - Bmag[:,yindex,:] = magnetic_field.Bmag(pol_grid.R, pol_grid.Z, ypos) - pressure[:,yindex,:] = magnetic_field.pressure(pol_grid.R, pol_grid.Z, ypos) + pol_grid, ypos = grid.getPoloidalGrid(yindex) + Bmag[:, yindex, :] = magnetic_field.Bmag(pol_grid.R, pol_grid.Z, ypos) + pressure[:, yindex, :] = magnetic_field.pressure(pol_grid.R, pol_grid.Z, ypos) + + metric["g_yy"][:, yindex, :] = (metric["g_yy"][:, yindex, :] + * (Bmag[:, yindex, :] + / magnetic_field.Byfunc(pol_grid.R, pol_grid.Z, ypos))**2) + metric["gyy"][:, yindex, :] = (metric["gyy"][:, yindex, :] + * (magnetic_field.Byfunc(pol_grid.R, pol_grid.Z, ypos) + / Bmag[:, yindex, :])**2) # Get attributes from magnetic field (e.g. psi) attributes = {} @@ -218,23 +221,24 @@ def write_maps(grid, magnetic_field, maps, gridfile='fci.grid.nc', attribute = np.zeros(grid.shape) for yindex in range(grid.numberOfPoloidalGrids()): pol_grid, ypos = grid.getPoloidalGrid(yindex) - attribute[:,yindex,:] = magnetic_field.attributes[name](pol_grid.R, pol_grid.Z, ypos) + attribute[:, yindex, :] = magnetic_field.attributes[name](pol_grid.R, pol_grid.Z, ypos) attributes[name] = attribute - + # Metric is now 3D if metric2d: # Remove the Z dimension from metric components - print("WARNING: Outputting 2D metrics, discarding metric information.") + if not quiet: + print("WARNING: Outputting 2D metrics, discarding metric information.") for key in metric: try: - metric[key] = metric[key][:,:,0] - except: + metric[key] = metric[key][:, :, 0] + except TypeError: pass # Make dz a constant - metric["dz"] = metric["dz"][0,0] + metric["dz"] = metric["dz"][0, 0] # Add Rxy, Bxy - metric["Rxy"] = maps["R"][:,:,0] - metric["Bxy"] = Bmag[:,:,0] + metric["Rxy"] = maps["R"][:, :, 0] + metric["Bxy"] = Bmag[:, :, 0] with bdata.DataFile(gridfile, write=True, create=True, format=format) as f: ixseps = nx+1 @@ -246,11 +250,11 @@ def write_maps(grid, magnetic_field, maps, gridfile='fci.grid.nc', f.write("dy", metric["dy"]) f.write("dz", metric["dz"]) - f.write("ixseps1",ixseps) - f.write("ixseps2",ixseps) + f.write("ixseps1", ixseps) + f.write("ixseps2", ixseps) # Metric tensor - + if new_names: for key, val in metric.items(): f.write(key, val) @@ -258,20 +262,20 @@ def write_maps(grid, magnetic_field, maps, gridfile='fci.grid.nc', # Translate between output variable names and metric names # Map from new to old names. Anything not in this dict # is output unchanged - name_changes = {"g_yy":"g_22", - "gyy":"g22", - "gxx":"g11", - "gxz":"g13", - "gzz":"g33", - "g_xx":"g_11", - "g_xz":"g_13", - "g_zz":"g_33"} + name_changes = {"g_yy": "g_22", + "gyy": "g22", + "gxx": "g11", + "gxz": "g13", + "gzz": "g33", + "g_xx": "g_11", + "g_xz": "g_13", + "g_zz": "g_33"} for key in metric: name = key if name in name_changes: name = name_changes[name] f.write(name, metric[key]) - + # Magnetic field f.write("B", Bmag) @@ -281,7 +285,7 @@ def write_maps(grid, magnetic_field, maps, gridfile='fci.grid.nc', # Attributes for name in attributes: f.write(name, attributes[name]) - + # Maps - write everything to file for key in maps: f.write(key, maps[key])