diff --git a/include/bout/paralleltransform.hxx b/include/bout/paralleltransform.hxx index aace91f21f..f61be64319 100644 --- a/include/bout/paralleltransform.hxx +++ b/include/bout/paralleltransform.hxx @@ -101,7 +101,13 @@ public: class ShiftedMetric : public ParallelTransform { public: ShiftedMetric() = delete; - ShiftedMetric(Mesh &mesh); + ShiftedMetric(Mesh &m) : mesh(m) {} + + /// initialize after constructor is finished, so that + /// localmesh->getCoordinateSystem() can call + /// ParallelTransform::getCoordinateSystem() for fields belonging to the + /// ParallelTransform + void initialize(); /*! * Calculates the yup() and ydown() fields of f @@ -151,6 +157,10 @@ private: Tensor yupPhs; ///< Cache of phase shifts for calculating yup fields Tensor ydownPhs; ///< Cache of phase shifts for calculating ydown fields +#if CHECK>0 + bool isinitialized = false; +#endif + /*! * Shift a 2D field in Z. * Since 2D fields are constant in Z, this has no effect diff --git a/src/mesh/mesh.cxx b/src/mesh/mesh.cxx index 2152052d30..80b9eaf64f 100644 --- a/src/mesh/mesh.cxx +++ b/src/mesh/mesh.cxx @@ -299,7 +299,8 @@ void Mesh::setParallelTransform() { }else if(ptstr == "shifted") { // Shifted metric method - transform = bout::utils::make_unique(*this); + transform = bout::utils::make_unique(*this); + static_cast(transform.get())->initialize(); }else if(ptstr == "fci") { @@ -307,7 +308,8 @@ void Mesh::setParallelTransform() { // Flux Coordinate Independent method bool fci_zperiodic; fci_options->get("z_periodic", fci_zperiodic, true); - transform = bout::utils::make_unique(*this, fci_zperiodic); + transform = bout::utils::make_unique(*this); + static_cast(transform.get())->initialize(fci_zperiodic); }else { throw BoutException(_("Unrecognised paralleltransform option.\n" diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx index 25002ffed3..6fe08dc9ac 100644 --- a/src/mesh/parallel/fci.cxx +++ b/src/mesh/parallel/fci.cxx @@ -53,8 +53,11 @@ inline BoutReal sgn(BoutReal val) { return (BoutReal(0) < val) - (val < BoutReal // Calculate all the coefficients needed for the spline interpolation // dir MUST be either +1 or -1 -FCIMap::FCIMap(Mesh &mesh, int dir, bool zperiodic) - : dir(dir), boundary_mask(mesh), corner_boundary_mask(mesh), y_prime(&mesh) { +void FCIMap::initialize(Mesh &mesh, int dir_in, bool zperiodic) { + dir = dir_in; + boundary_mask = BoutMask(mesh); + corner_boundary_mask = BoutMask(mesh); + y_prime = Field3D(0., &mesh); interp = InterpolationFactory::getInstance()->create(&mesh); interp->setYOffset(dir); @@ -238,7 +241,7 @@ FCIMap::FCIMap(Mesh &mesh, int dir, bool zperiodic) const Field3D FCIMap::integrate(Field3D &f) const { TRACE("FCIMap::integrate"); - + // Cell centre values Field3D centre = interp->interpolate(f); @@ -296,6 +299,8 @@ const Field3D FCIMap::integrate(Field3D &f) const { void FCITransform::calcYUpDown(Field3D &f) { TRACE("FCITransform::calcYUpDown"); + ASSERT1(isinitialized); + // Ensure that yup and ydown are different fields f.splitYupYdown(); @@ -306,7 +311,9 @@ void FCITransform::calcYUpDown(Field3D &f) { void FCITransform::integrateYUpDown(Field3D &f) { TRACE("FCITransform::integrateYUpDown"); - + + ASSERT1(isinitialized); + // Ensure that yup and ydown are different fields f.splitYupYdown(); diff --git a/src/mesh/parallel/fci.hxx b/src/mesh/parallel/fci.hxx index 6e3e7c90a3..4eab67c13e 100644 --- a/src/mesh/parallel/fci.hxx +++ b/src/mesh/parallel/fci.hxx @@ -39,12 +39,10 @@ class FCIMap { /// Interpolation object Interpolation *interp; // Cell centre Interpolation *interp_corner; // Cell corner at (x+1, z+1) - - /// Private constructor - must be initialised with mesh - FCIMap(); public: + /// initialize method must be called before using FCIMap /// dir MUST be either +1 or -1 - FCIMap(Mesh& mesh, int dir, bool zperiodic); + void initialize(Mesh& mesh, int dir_in, bool zperiodic_in); int dir; /**< Direction of map */ @@ -65,9 +63,21 @@ public: */ 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(Mesh& m) : mesh(m) {} + + /// initialize after constructor is finished, so that + /// localmesh->getCoordinateSystem() can call + /// ParallelTransform::getCoordinateSystem() for fields belonging to the + /// ParallelTransform + void initialize(bool zperiodic_in = true) { + forward_map.initialize(mesh, +1, zperiodic); + backward_map.initialize(mesh, -1, zperiodic); + zperiodic = zperiodic_in; + +#if CHECK>0 + isinitialized = true; +#endif + } void calcYUpDown(Field3D &f) override; @@ -97,6 +107,10 @@ private: FCIMap backward_map; /**< FCI map for field lines in -ve y */ bool zperiodic; /**< Is the z-direction periodic? */ + +#if CHECK>0 + bool isinitialized = false; +#endif }; #endif // __FCITRANSFORM_H__ diff --git a/src/mesh/parallel/shiftedmetric.cxx b/src/mesh/parallel/shiftedmetric.cxx index b3ddcb8676..1f960372a5 100644 --- a/src/mesh/parallel/shiftedmetric.cxx +++ b/src/mesh/parallel/shiftedmetric.cxx @@ -15,9 +15,9 @@ #include -ShiftedMetric::ShiftedMetric(Mesh &m) : mesh(m), zShift(&m) { +void ShiftedMetric::initialize() { // Read the zShift angle from the mesh - + zShift = Field2D(&mesh); if(mesh.get(zShift, "zShift")) { // No zShift variable. Try qinty in BOUT grid files mesh.get(zShift, "qinty"); @@ -76,6 +76,9 @@ ShiftedMetric::ShiftedMetric(Mesh &m) : mesh(m), zShift(&m) { } } +#if CHECK > 1 + isinitialized = true; +#endif } /*! @@ -83,6 +86,7 @@ ShiftedMetric::ShiftedMetric(Mesh &m) : mesh(m), zShift(&m) { */ void ShiftedMetric::calcYUpDown(Field3D &f) { ASSERT1(&mesh == f.getMesh()); + ASSERT1(isinitialized); f.splitYupYdown(); @@ -106,7 +110,9 @@ void ShiftedMetric::calcYUpDown(Field3D &f) { * and Y is then field aligned. */ const Field3D ShiftedMetric::toFieldAligned(const Field3D &f, const REGION region) { + ASSERT1(isinitialized); ASSERT2(f.getCoordinateSystem() == COORDINATE_SYSTEM::Orthogonal); + Field3D result = shiftZ(f, toAlignedPhs, region); result.setCoordinateSystem(COORDINATE_SYSTEM::FieldAligned); return result; @@ -117,7 +123,9 @@ const Field3D ShiftedMetric::toFieldAligned(const Field3D &f, const REGION regio * but Y is not field aligned. */ const Field3D ShiftedMetric::fromFieldAligned(const Field3D &f, const REGION region) { + ASSERT1(isinitialized); ASSERT2(f.getCoordinateSystem() == COORDINATE_SYSTEM::FieldAligned); + Field3D result = shiftZ(f, fromAlignedPhs, region); result.setCoordinateSystem(f.getMesh()->getCoordinateSystem()); return result; @@ -126,6 +134,7 @@ const Field3D ShiftedMetric::fromFieldAligned(const Field3D &f, const REGION reg const Field3D ShiftedMetric::shiftZ(const Field3D& f, const Tensor& phs, const REGION region) { ASSERT1(&mesh == f.getMesh()); + ASSERT1(isinitialized); if(mesh.LocalNz == 1) return f; // Shifting makes no difference @@ -141,6 +150,8 @@ const Field3D ShiftedMetric::shiftZ(const Field3D& f, const Tensor& ph } void ShiftedMetric::shiftZ(const BoutReal* in, const dcomplex* phs, BoutReal* out) { + ASSERT1(isinitialized); + Array cmplx(nmodes); // Take forward FFT @@ -160,6 +171,7 @@ void ShiftedMetric::shiftZ(const BoutReal* in, const dcomplex* phs, BoutReal* ou //Old approach retained so we can still specify a general zShift const Field3D ShiftedMetric::shiftZ(const Field3D &f, const Field2D &zangle, const REGION region) { + ASSERT1(isinitialized); ASSERT1(&mesh == f.getMesh()); ASSERT1(f.getLocation() == zangle.getLocation()); if(mesh.LocalNz == 1) @@ -182,6 +194,8 @@ const Field3D ShiftedMetric::shiftZ(const Field3D &f, const Field2D &zangle, con } void ShiftedMetric::shiftZ(const BoutReal *in, int len, BoutReal zangle, BoutReal *out) { + ASSERT1(isinitialized); + int nmodes = len/2 + 1; // Complex array used for FFTs diff --git a/tests/integrated/test-yupdown/test_yupdown.cxx b/tests/integrated/test-yupdown/test_yupdown.cxx index cd5ea6ed3b..604d31aa44 100644 --- a/tests/integrated/test-yupdown/test_yupdown.cxx +++ b/tests/integrated/test-yupdown/test_yupdown.cxx @@ -39,6 +39,7 @@ int main(int argc, char** argv) { BoutInitialise(argc, argv); ShiftedMetric s(*mesh); + s.initialize(); // Read variable from mesh Field3D var;