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

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.

int invert_laplace(const FieldPerp &b, FieldPerp &x, int flags, const Field2D *a, const Field2D *c = nullptr, const Field2D *d = nullptr)
int invert_laplace(const Field3D &b, Field3D &x, int flags, const Field2D *a, const Field2D *c = nullptr, const Field2D *d = nullptr)
const Field3D invert_laplace(const Field3D &b, int flags, const Field2D *a = nullptr, const Field2D *c = nullptr, const Field2D *d = nullptr)

More readable API for calling Laplacian inversion. Returns x.

Variables

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 Laplacian
#include <invert_laplace.hxx>

Base class for Laplacian inversion.

Subclassed by LaplaceCyclic, LaplaceMultigrid, LaplaceMumps, LaplaceNaulin, LaplacePDD, LaplacePetsc, LaplaceSerialBand, LaplaceSerialTri, LaplaceShoot, LaplaceSPT

Public Functions

Laplacian(Options *options = nullptr, const CELL_LOC loc = CELL_CENTRE, Mesh *mesh_in = nullptr)

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

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

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

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

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

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

Performs the laplacian inversion y-slice by y-slice

Return
x All the y-slices of x_slice in the equation A*x_slice = b_slice
Parameters
  • b: All the y-slices of b_slice, which is the right hand side of the equation A*x_slice = b_slice
  • x0: All the y-slices of the variable eventually used to set BC

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

Laplacian *create(Options *opt = nullptr, const CELL_LOC loc = CELL_CENTRE, Mesh *mesh_in = nullptr)

Create a new Laplacian solver

Parameters
  • opt: The options section to use. By default “laplace” will be used

Laplacian *defaultInstance()

Return pointer to global singleton.

void cleanup()

Frees all memory.

Protected Functions

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)
void tridagMatrix(dcomplex **avec, dcomplex **bvec, dcomplex **cvec, dcomplex **bk, int jy, int flags, int inner_boundary_flags, int outer_boundary_flags, const Field2D *a = nullptr, const Field2D *ccoef = nullptr, const Field2D *d = nullptr)

Sets the coefficients for parallel tridiagonal matrix inversion.

Uses the laplace_tridag_coefs routine above to fill a matrix [kz][ix] of coefficients

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: Lower diagonal of the tridiagonal matrix. DO NOT CONFUSE WITH “A”
  • bvec: The main diagonal
  • cvec: The upper diagonal. DO NOT CONFUSE WITH “C” (called ccoef here)
  • bk: The b in Ax = b
  • jy: Index of the current y-slice
  • kz: The mode number index
  • kwave: 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: Global flags of the inversion
  • inner_boundary_flags: Flags used to set the inner boundary
  • outer_boundary_flags: Flags used to set the outer boundary
  • a: A in the equation above. DO NOT CONFUSE WITH avec
  • c1coef: C1 in the equation above. DO NOT CONFUSE WITH cvec
  • c2coef: C2 in the equation above. DO NOT CONFUSE WITH cvec
  • d: D in the equation above
  • includeguards: Whether or not the guard points in x should be used
  • avec: Lower diagonal of the tridiagonal matrix. DO NOT CONFUSE WITH “A”
  • bvec: The main diagonal
  • cvec: 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

Laplacian *instance = nullptr

Singleton instance.