File invert_laplace.hxx

Perpendicular Laplacian inversion using FFT and Tridiagonal solver

Equation solved is: \(d*\nabla^2_\perp x + (1/c)\nabla_perp c\cdot\nabla_\perp x + a x = b\)

Where a, c and d are functions of x and y only (not z)

Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu

Contact: Ben Dudson, bd512@york.ac.uk

This file is part of BOUT++.

BOUT++ is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

BOUT++ is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License along with BOUT++. If not, see http://www.gnu.org/licenses/.

Defines

PVEC_REAL_MPI_TYPE

Typedefs

using RegisterLaplace = LaplaceFactory::RegisterInFactory<DerivedType>

Simpler name for Factory registration helper class

Usage:

#include <bout/laplacefactory.hxx>
namespace {
RegisterLaplace<MyLaplace> registerlaplacemine("mylaplace");
}

using RegisterUnavailableLaplace = LaplaceFactory::RegisterUnavailableInFactory

Functions

void laplace_tridag_coefs(int jx, int jy, int jz, dcomplex &a, dcomplex &b, dcomplex &c, const Field2D *ccoef = nullptr, const Field2D *d = nullptr, CELL_LOC loc = CELL_DEFAULT)

Returns the coefficients for a tridiagonal matrix for laplace. Used by Delp2 too.

Variables

constexpr auto LAPLACE_SPT = "spt"
constexpr auto LAPLACE_PDD = "pdd"
constexpr auto LAPLACE_TRI = "tri"
constexpr auto LAPLACE_BAND = "band"
constexpr auto LAPLACE_PETSC = "petsc"
constexpr auto LAPLACE_PETSCAMG = "petscamg"
constexpr auto LAPLACE_PETSC3DAMG = "petsc3damg"
constexpr auto LAPLACE_HYPRE3D = "hypre3d"
constexpr auto LAPLACE_CYCLIC = "cyclic"
constexpr auto LAPLACE_MULTIGRID = "multigrid"
constexpr auto LAPLACE_NAULIN = "naulin"
constexpr auto LAPLACE_IPT = "ipt"
constexpr auto LAPLACE_PCR = "pcr"
constexpr int INVERT_DC_GRAD = 1

Zero-gradient for DC (constant in Z) component. Default is zero value.

constexpr int INVERT_AC_GRAD = 2

Zero-gradient for AC (non-constant in Z) component. Default is zero value.

constexpr int INVERT_AC_LAP = 4

Use zero-laplacian (decaying solution) to AC component.

constexpr int INVERT_SYM = 8

Use symmetry to enforce either zero-value or zero-gradient.

constexpr int INVERT_SET = 16

Set boundary to value.

constexpr int INVERT_RHS = 32

Use input value in RHS boundary.

constexpr int INVERT_DC_LAP = 64

Use zero-laplacian solution for DC component.

constexpr int INVERT_BNDRY_ONE = 128

Only use one boundary point.

constexpr int INVERT_DC_GRADPAR = 256
constexpr int INVERT_DC_GRADPARINV = 512
constexpr int INVERT_IN_CYLINDER = 1024

For use in cylindrical coordiate system.

constexpr int INVERT_ZERO_DC = 1

Zero the DC (constant in Z) component of the solution.

constexpr int INVERT_START_NEW = 2

Iterative method start from solution=0. Has no effect for direct solvers.

constexpr int INVERT_BOTH_BNDRY_ONE = 4

Sets the width of the boundaries to 1.

constexpr int INVERT_4TH_ORDER = 8

Use band solver for 4th order in x.

constexpr int INVERT_KX_ZERO = 16

Zero the kx=0, n = 0 component.

class LaplaceFactory : public Factory<Laplacian, LaplaceFactory, Options*, CELL_LOC, Mesh*, Solver*, Datafile*>

Public Functions

inline ReturnType create(Options *options = nullptr, CELL_LOC loc = CELL_CENTRE, Mesh *mesh = nullptr, Solver *solver = nullptr, Datafile *dump = nullptr)
inline ReturnType create(const std::string &type, Options *options) const

Public Static Attributes

static constexpr auto type_name = "Laplacian"
static constexpr auto section_name = "laplace"
static constexpr auto option_name = "type"
static constexpr auto default_type = LAPLACE_CYCLIC
class Laplacian
#include <invert_laplace.hxx>

Base class for Laplacian inversion.

Subclassed by LaplaceCyclic, LaplaceIPT, LaplaceMultigrid, LaplaceNaulin, LaplacePCR, LaplacePDD, LaplacePetsc, LaplacePetsc3dAmg, LaplaceSerialBand, LaplaceSerialTri, LaplaceSPT

Public Functions

Laplacian(Options *options = nullptr, const CELL_LOC loc = CELL_CENTRE, Mesh *mesh_in = nullptr, Solver *solver = nullptr, Datafile *dump = nullptr)

Laplacian inversion initialisation. Called once at the start to get settings.

virtual ~Laplacian() = default
virtual void setCoefA(const Field2D &val) = 0

Set coefficients for inversion. Re-builds matrices if necessary.

inline virtual void setCoefA(const Field3D &val)
inline virtual void setCoefA(BoutReal r)
virtual void setCoefC(const Field2D &val) = 0
inline virtual void setCoefC(const Field3D &val)
inline virtual void setCoefC(BoutReal r)
inline virtual void setCoefC1(const Field2D &val)
inline virtual void setCoefC1(const Field3D &val)
inline virtual void setCoefC1(BoutReal r)
inline virtual void setCoefC2(const Field2D &val)
inline virtual void setCoefC2(const Field3D &val)
inline virtual void setCoefC2(BoutReal r)
virtual void setCoefD(const Field2D &val) = 0
inline virtual void setCoefD(const Field3D &val)
inline virtual void setCoefD(BoutReal r)
virtual void setCoefEx(const Field2D &val) = 0
inline virtual void setCoefEx(const Field3D &val)
inline virtual void setCoefEx(BoutReal r)
virtual void setCoefEz(const Field2D &val) = 0
inline virtual void setCoefEz(const Field3D &val)
inline virtual void setCoefEz(BoutReal r)
inline virtual void setGlobalFlags(int f)
inline virtual void setInnerBoundaryFlags(int f)
inline virtual void setOuterBoundaryFlags(int f)
inline virtual bool uses3DCoefs() const

Does this solver use Field3D coefficients (true) or only their DC component (false)

virtual FieldPerp solve(const FieldPerp &b) = 0
virtual Field3D solve(const Field3D &b)
virtual Field2D solve(const Field2D &b)
inline virtual FieldPerp solve(const FieldPerp &b, const FieldPerp &x0)
virtual Field3D solve(const Field3D &b, const Field3D &x0)

Performs the laplacian inversion y-slice by y-slice

Parameters
  • b[in] All the y-slices of b_slice, which is the right hand side of the equation A*x_slice = b_slice

  • x0[in] All the y-slices of the variable eventually used to set BC

Returns

x All the y-slices of x_slice in the equation A*x_slice = b_slice

virtual Field2D solve(const Field2D &b, const Field2D &x0)
void tridagCoefs(int jx, int jy, int jz, dcomplex &a, dcomplex &b, dcomplex &c, const Field2D *ccoef = nullptr, const Field2D *d = nullptr, CELL_LOC loc = CELL_DEFAULT)

Coefficients in tridiagonal inversion.

Public Static Functions

static inline std::unique_ptr<Laplacian> create(Options *opts = nullptr, const CELL_LOC location = CELL_CENTRE, Mesh *mesh_in = nullptr, Solver *solver = nullptr, Datafile *dump = nullptr)

Create a new Laplacian solver

Parameters

opt[in] The options section to use. By default “laplace” will be used

static Laplacian *defaultInstance()

Return pointer to global singleton.

static void cleanup()

Frees all memory.

Protected Functions

inline void tridagCoefs(int jx, int jy, BoutReal kwave, dcomplex &a, dcomplex &b, dcomplex &c, const Field2D *ccoef = nullptr, const Field2D *d = nullptr, CELL_LOC loc = CELL_DEFAULT)
void tridagCoefs(int jx, int jy, BoutReal kwave, dcomplex &a, dcomplex &b, dcomplex &c, const Field2D *c1coef, const Field2D *c2coef, const Field2D *d, CELL_LOC loc = CELL_DEFAULT)
inline void tridagMatrix(dcomplex *avec, dcomplex *bvec, dcomplex *cvec, dcomplex *bk, int jy, int kz, BoutReal kwave, int flags, int inner_boundary_flags, int outer_boundary_flags, const Field2D *a, const Field2D *ccoef, const Field2D *d, bool includeguards = true)
void tridagMatrix(dcomplex *avec, dcomplex *bvec, dcomplex *cvec, dcomplex *bk, int jy, int kz, BoutReal kwave, int flags, int inner_boundary_flags, int outer_boundary_flags, const Field2D *a, const Field2D *c1coef, const Field2D *c2coef, const Field2D *d, bool includeguards = true)

Set the matrix components of A in Ax=b

This function will

  1. Calling tridagCoef, solving

    D*Laplace_perp(x) + (1/C1)Grad_perp(C2)*Grad_perp(x) + Ax = B

    for each fourier component

  2. Set the boundary conditions by setting the first and last rows properly

Parameters
  • avec[in] Lower diagonal of the tridiagonal matrix. DO NOT CONFUSE WITH “A”

  • bvec[in] The main diagonal

  • cvec[in] The upper diagonal. DO NOT CONFUSE WITH “C” (called ccoef here)

  • bk[in] The b in Ax = b

  • jy[in] Index of the current y-slice

  • kz[in] The mode number index

  • kwave[in] The mode number (different from kz only if we are taking a part of the z-domain [and not from 0 to 2*pi])

  • global_flags[in] Global flags of the inversion

  • inner_boundary_flags[in] Flags used to set the inner boundary

  • outer_boundary_flags[in] Flags used to set the outer boundary

  • a[in] A in the equation above. DO NOT CONFUSE WITH avec

  • c1coef[in] C1 in the equation above. DO NOT CONFUSE WITH cvec

  • c2coef[in] C2 in the equation above. DO NOT CONFUSE WITH cvec

  • d[in] D in the equation above

  • includeguards[in] Whether or not the guard points in x should be used

  • avec[out] Lower diagonal of the tridiagonal matrix. DO NOT CONFUSE WITH “A”

  • bvec[out] The main diagonal

  • cvec[out] The upper diagonal. DO NOT CONFUSE WITH “C” (called ccoef here)

Protected Attributes

bool async_send

If true, use asyncronous send in parallel algorithms.

int maxmode

The maximum Z mode to solve for.

bool low_mem

If true, reduce the amount of memory used.

bool all_terms

applies to Delp2 operator and laplacian inversion

bool nonuniform

Non-uniform mesh correction.

bool include_yguards

solve in y-guard cells, default true.

int extra_yguards_lower

exclude some number of points at the lower boundary, useful for staggered grids or when boundary conditions make inversion redundant

int extra_yguards_upper

exclude some number of points at the upper boundary, useful for staggered grids or when boundary conditions make inversion redundant

int global_flags

Default flags.

int inner_boundary_flags

Flags to set inner boundary condition.

int outer_boundary_flags

Flags to set outer boundary condition.

CELL_LOC location

staggered grid location of this solver

Mesh *localmesh

Mesh object for this solver.

Coordinates *coords

Coordinates object, so we only have to call localmesh->getCoordinates(location) once

Private Static Attributes

static std::unique_ptr<Laplacian> instance = nullptr

Singleton instance.