dg::Elliptic3d< Geometry, Matrix, Container > Class Template Reference

A 3d negative elliptic differential operator \( -\nabla \cdot ( \mathbf{\chi}\cdot \nabla ) \). More...

Public Types
using	geometry_type = Geometry
using	matrix_type = Matrix
using	container_type = Container
using	value_type = get_value_type< Container >

Public Member Functions
	Elliptic3d ()=default
	empty object ( no memory allocation) More...
	Elliptic3d (const Geometry &g, direction dir=forward, value_type jfactor=1., bool chi_weight_jump=false)
	Construct from Grid. More...
	Elliptic3d (const Geometry &g, bc bcx, bc bcy, bc bcz, direction dir=forward, value_type jfactor=1., bool chi_weight_jump=false)
	Construct from grid and boundary conditions. More...
template<class ... Params>
void	construct (Params &&...ps)
	Perfect forward parameters to one of the constructors. More...
template<class ContainerType0 >
void	set_chi (const ContainerType0 &sigma)
	Change scalar part in Chi tensor. More...
template<class ContainerType0 >
void	set_chi (const SparseTensor< ContainerType0 > &tau)
	Change tensor part in Chi tensor. More...
const Container &	weights () const
	Return the weights making the operator self-adjoint. More...
const Container &	precond () const
	Return the default preconditioner to use in conjugate gradient. More...
void	set_jfactor (value_type new_jfactor)
	Set the currently used jfactor ( \( \alpha \)) More...
value_type	get_jfactor () const
	Get the currently used jfactor ( \( \alpha \)) More...
void	set_jump_weighting (bool jump_weighting)
	Set the chi weighting of jump terms. More...
bool	get_jump_weighting () const
	Get the current state of chi weighted jump terms. More...
void	set_compute_in_2d (bool compute_in_2d)
	Restrict the problem to the first 2 dimensions. More...
template<class ContainerType0 , class ContainerType1 >
void	symv (const ContainerType0 &x, ContainerType1 &y)
	Compute elliptic term and store in output. More...
template<class ContainerType0 , class ContainerType1 >
void	symv (value_type alpha, const ContainerType0 &x, value_type beta, ContainerType1 &y)
	Compute elliptic term and add to output. More...
template<class ContainerType0 , class ContainerType1 >
void	variation (const ContainerType0 &phi, ContainerType1 &sigma)
	\( \sigma = (\nabla\phi\cdot\tau\cdot\nabla \phi) \) More...
template<class ContainerTypeL , class ContainerType0 , class ContainerType1 >
void	variation (const ContainerTypeL &lambda, const ContainerType0 &phi, ContainerType1 &sigma)
	\( \sigma = \lambda^2(\nabla\phi\cdot\tau\cdot\nabla \phi) \) More...
template<class ContainerTypeL , class ContainerType0 , class ContainerType1 >
void	variation (value_type alpha, const ContainerTypeL &lambda, const ContainerType0 &phi, value_type beta, ContainerType1 &sigma)
	\( \sigma = \alpha \lambda^2 (\nabla\phi\cdot\tau\cdot\nabla \phi) + \beta \sigma\) More...

Detailed Description

template<class Geometry, class Matrix, class Container>
class dg::Elliptic3d< Geometry, Matrix, Container >

A 3d negative elliptic differential operator \( -\nabla \cdot ( \mathbf{\chi}\cdot \nabla ) \).

The term discretized is

\[ -\nabla \cdot ( \mathbf{\chi} \cdot \nabla ) \]

where \( \nabla \) is the two-dimensional nabla and \(\chi = \sigma \mathbf{\tau}\) is a (possibly spatially dependent) tensor with scalar part \( \sigma\) (usually the volume form) and tensor part \( \tau\) (usually the inverse metric). In general coordinates that means

\[ -\frac{1}{\sqrt{g}}\left( \partial_x\left(\sqrt{g}\left(\chi^{xx}\partial_x + \chi^{xy}\partial_y + \chi^{xz}\partial_z \right)\right) + \partial_y\left(\sqrt{g}\left(\chi^{yx}\partial_x + \chi^{yy}\partial_y + \chi^{yz}\partial_z \right)\right) + \partial_z\left(\sqrt{g}\left(\chi^{zx}\partial_x + \chi^{zy}\partial_y + \chi^{zz}\partial_z \right)\right) \right)\]

is discretized. Per default, \( \chi = \sqrt{g} g^{-1}\) but you can set it to any tensor you like (in order for the operator to be invertible \(\chi\) should be symmetric and positive definite though).

See also

Our theory guide Introduction to dg methods on overleaf holds a detailed derivation

Note that the local discontinuous Galerkin discretization adds so-called jump terms

\[ D^\dagger \chi D + \alpha\chi_{on/off} J \]

where \(\alpha\) is a scale factor ( = jfactor), \( D \) contains the discretizations of the above derivatives, and \( J\) is a self-adjoint matrix. ( \(J\) is added before the volume element is divided). The adjoint of a matrix is defined with respect to the volume element including dG weights. Usually the default \( \alpha=1 \) is a good choice. However, in some cases, e.g. when \( \sigma \) exhibits very large variations \( \alpha=0.1\) or \( \alpha=0.01\) might be better values. In a time dependent problem the value of \(\alpha\) determines the numerical diffusion, i.e. for too low values numerical oscillations may appear. The \( \chi_{on/off} \) in the jump term serves to weight the jump term with \( \chi \). This can be switched either on or off with off being the default. Also note that a forward discretization has more diffusion than a centered discretization.

The following code snippet demonstrates the use of Elliptic3d in an inversion problem

    dg::CylindricalGrid3d grid( R_0, R_0+lx, 0, ly, 0,lz, n, Nx, Ny,Nz, bcx, bcy, bcz);
    dg::DVec w3d = dg::create::volume( grid);
    dg::DVec x = dg::evaluate( initial, grid);
 
    dg::Elliptic3d<dg::aGeometry3d, dg::DMatrix, dg::DVec> laplace(grid, dg::centered);
 
    laplace.set_jump_weighting(jump_weight);
 
    dg::PCG< dg::DVec > pcg( x, n*n*Nx*Ny*Nz);
 
    const dg::DVec solution = dg::evaluate ( fct, grid);
    dg::DVec b = dg::evaluate ( laplace3d_fct, grid);
 
    std::cout << "For a precision of "<< eps<<" ..."<<std::endl;
    unsigned num;
    t.tic();
    num = pcg.solve( laplace, x, b, 1., w3d, eps);
    t.toc();
    std::cout << "Number of pcg iterations "<<num<<std::endl;
    std::cout << "... on the device took   "<< t.diff()<<"s\n";

Template Parameters

Geometry	A type that is or derives from one of the abstract geometry base classes ( aGeometry2d, aGeometry3d, aMPIGeometry2d, ...). Geometry determines which Matrix and Container types can be used:
Matrix	A class for which the dg::blas2::symv functions are callable in connection with the Container class and to which the return type of dg::create::dx() can be converted using dg::blas2::transfer. The Matrix type can be one of: dg::HMatrix with dg::HVec and one of the shared memory geometries dg::DMatrix with dg::DVec and one of the shared memory geometries dg::MHMatrix with dg::MHVec and one of the MPI geometries dg::MDMatrix with dg::MDVec and one of the MPI geometries
Container	A data container class for which the blas1 functionality is overloaded and to which the return type of blas1::subroutine() can be converted using dg::assign. We assume that Container is copyable/assignable and has a swap member function. In connection with Geometry this is one of dg::HVec, dg::DVec when Geometry is a shared memory geometry dg::MHVec or dg::MDVec when Geometry is one of the MPI geometries This class has the SelfMadeMatrixTag so it can be used in blas2::symv functions and thus in a conjugate gradient solver.

Note

The constructors initialize \( \chi=\sqrt{g}g^{-1}\) so that a negative laplacian operator results

The inverse of \( \sigma\) makes a good general purpose preconditioner

Since the pattern arises quite often (because of the ExB velocity \( u_E^2\) in the ion gyro-centre potential) this class also can compute the variation integrand \( \lambda^2\nabla \phi\cdot \tau\cdot\nabla\phi\)

Attention

Pay attention to the negative sign which is necessary to make the matrix positive definite

Member Typedef Documentation

container_type

template<class Geometry , class Matrix , class Container >

using dg::Elliptic3d< Geometry, Matrix, Container >::container_type = Container

geometry_type

template<class Geometry , class Matrix , class Container >

using dg::Elliptic3d< Geometry, Matrix, Container >::geometry_type = Geometry

matrix_type

template<class Geometry , class Matrix , class Container >

using dg::Elliptic3d< Geometry, Matrix, Container >::matrix_type = Matrix

value_type

template<class Geometry , class Matrix , class Container >

using dg::Elliptic3d< Geometry, Matrix, Container >::value_type = get_value_type<Container>

Constructor & Destructor Documentation

Elliptic3d() [1/3]

template<class Geometry , class Matrix , class Container >

dg::Elliptic3d< Geometry, Matrix, Container >::Elliptic3d ( )

default

empty object ( no memory allocation)

Elliptic3d() [2/3]

template<class Geometry , class Matrix , class Container >

dg::Elliptic3d< Geometry, Matrix, Container >::Elliptic3d ( const Geometry &  g, direction  dir = forward, value_type  jfactor = 1., bool  chi_weight_jump = false  )

inline

Construct from Grid.

Parameters

g	The Grid; boundary conditions are taken from here
dir	Direction of the right first derivative in x and y (i.e. dg::forward, dg::backward or dg::centered), the direction of the z derivative is dg::centered if nz = 1
jfactor	( \( = \alpha \) ) scale jump terms (1 is a good value but in some cases 0.1 or 0.01 might be better)
chi_weight_jump	If true, the Jump terms are multiplied with the Chi matrix, else it is ignored

Note

chi is assumed the metric per default

Elliptic3d() [3/3]

template<class Geometry , class Matrix , class Container >

dg::Elliptic3d< Geometry, Matrix, Container >::Elliptic3d ( const Geometry &  g, bc  bcx, bc  bcy, bc  bcz, direction  dir = forward, value_type  jfactor = 1., bool  chi_weight_jump = false  )

inline

Construct from grid and boundary conditions.

Parameters

g	The Grid
bcx	boundary condition in x
bcy	boundary contition in y
bcz	boundary contition in z
dir	Direction of the right first derivative in x and y (i.e. dg::forward, dg::backward or dg::centered), the direction of the z derivative is always dg::centered
jfactor	( \( = \alpha \) ) scale jump terms (1 is a good value but in some cases 0.1 or 0.01 might be better)
chi_weight_jump	If true, the Jump terms are multiplied with the Chi matrix, else it is ignored

Note

chi is the metric tensor multiplied by the volume element per default

Member Function Documentation

construct()

template<class Geometry , class Matrix , class Container >

template<class ... Params>

void dg::Elliptic3d< Geometry, Matrix, Container >::construct ( Params &&...  ps )

inline

Perfect forward parameters to one of the constructors.

Template Parameters

Params

deduced by the compiler

Parameters

ps	parameters forwarded to constructors

get_jfactor()

template<class Geometry , class Matrix , class Container >

value_type dg::Elliptic3d< Geometry, Matrix, Container >::get_jfactor ( ) const

inline

Get the currently used jfactor ( \( \alpha \))

Returns

The current scale factor for jump terms

get_jump_weighting()

template<class Geometry , class Matrix , class Container >

bool dg::Elliptic3d< Geometry, Matrix, Container >::get_jump_weighting ( ) const

inline

Get the current state of chi weighted jump terms.

Returns

Whether the weighting of jump terms with chi is enabled. Either true or false.

precond()

template<class Geometry , class Matrix , class Container >

const Container & dg::Elliptic3d< Geometry, Matrix, Container >::precond ( ) const

inline

Return the default preconditioner to use in conjugate gradient.

Currently returns 1 divided by the scalar part of \( \sigma\). This is especially good when \( \sigma\) exhibits large amplitudes or variations

Returns

the inverse of \(\sigma\).

set_chi() [1/2]

template<class Geometry , class Matrix , class Container >

template<class ContainerType0 >

void dg::Elliptic3d< Geometry, Matrix, Container >::set_chi ( const ContainerType0 &  sigma )

inline

Change scalar part in Chi tensor.

Internally, we split the tensor \(\chi = \sigma\mathbf{\tau}\) into a scalar part \( \sigma\) and a tensor part \( \tau\) and you can set each part seperately. This functions sets the scalar part.

Parameters

sigma

The new scalar part in \(\chi\)

Note

The class will take care of the volume element in the divergence so do not multiply it to sigma yourself

Attention

If some or all elements of sigma are zero the preconditioner is invalidated and the operator can no longer be inverted (due to divide by zero) until set_chi is called with a positive sigma again. The symv function can always be called, however, if sigma is zero you likely also want to set the jfactor to 0 because the jumps in phi may not vanish and then pollute the result.

Template Parameters

ContainerType0

must be usable with Container in The dg dispatch system

set_chi() [2/2]

template<class Geometry , class Matrix , class Container >

template<class ContainerType0 >

void dg::Elliptic3d< Geometry, Matrix, Container >::set_chi ( const SparseTensor< ContainerType0 > &  tau )

inline

Change tensor part in Chi tensor.

We split the tensor \(\chi = \sigma\mathbf{\tau}\) into a scalar part \( \sigma\) and a tensor part \( \tau\) and you can set each part seperately. This functions sets the tensor part.

Note

The class will take care of the volume element in the divergence so do not multiply it to tau yourself

Parameters

tau	The new tensor part in \(\chi\) (must be positive definite)

Note

the 3d parts in tau will be ignored for 2d computations

Template Parameters

ContainerType0

must be usable in dg::assign to Container

set_compute_in_2d()

template<class Geometry , class Matrix , class Container >

void dg::Elliptic3d< Geometry, Matrix, Container >::set_compute_in_2d ( bool  compute_in_2d )

inline

Restrict the problem to the first 2 dimensions.

This effectively makes the behaviour of dg::Elliptic3d identical to the dg::Elliptic class.

Parameters

compute_in_2d

if true, the symv function avoids the derivative in z, false reverts to the original behaviour.

set_jfactor()

template<class Geometry , class Matrix , class Container >

void dg::Elliptic3d< Geometry, Matrix, Container >::set_jfactor ( value_type  new_jfactor )

inline

Set the currently used jfactor ( \( \alpha \))

Parameters

new_jfactor

The new scale factor for jump terms

set_jump_weighting()

template<class Geometry , class Matrix , class Container >

void dg::Elliptic3d< Geometry, Matrix, Container >::set_jump_weighting ( bool  jump_weighting )

inline

Set the chi weighting of jump terms.

Parameters

jump_weighting

Switch for weighting the jump factor with chi. Either true or false.

symv() [1/2]

template<class Geometry , class Matrix , class Container >

template<class ContainerType0 , class ContainerType1 >

void dg::Elliptic3d< Geometry, Matrix, Container >::symv ( const ContainerType0 &  x, ContainerType1 &  y  )

inline

Compute elliptic term and store in output.

i.e. y=M*x

Parameters

x	left-hand-side
y	result

Template Parameters

ContainerTypes

must be usable with Container in The dg dispatch system

symv() [2/2]

template<class Geometry , class Matrix , class Container >

template<class ContainerType0 , class ContainerType1 >

void dg::Elliptic3d< Geometry, Matrix, Container >::symv ( value_type  alpha, const ContainerType0 &  x, value_type  beta, ContainerType1 &  y  )

inline

Compute elliptic term and add to output.

i.e. y=alpha*M*x+beta*y

Parameters

alpha	a scalar
x	left-hand-side
beta	a scalar
y	result

Template Parameters

ContainerTypes

must be usable with Container in The dg dispatch system

variation() [1/3]

template<class Geometry , class Matrix , class Container >

template<class ContainerType0 , class ContainerType1 >

void dg::Elliptic3d< Geometry, Matrix, Container >::variation ( const ContainerType0 &  phi, ContainerType1 &  sigma  )

inline

\( \sigma = (\nabla\phi\cdot\tau\cdot\nabla \phi) \)

where \( \tau \) is the tensor part of \(\chi\) that is the (inverse) metric by default

Parameters

phi	the vector to take the variation of
sigma	(inout) the variation

Template Parameters

ContainerTypes

must be usable with Container in The dg dispatch system

variation() [2/3]

template<class Geometry , class Matrix , class Container >

template<class ContainerTypeL , class ContainerType0 , class ContainerType1 >

void dg::Elliptic3d< Geometry, Matrix, Container >::variation ( const ContainerTypeL &  lambda, const ContainerType0 &  phi, ContainerType1 &  sigma  )

inline

\( \sigma = \lambda^2(\nabla\phi\cdot\tau\cdot\nabla \phi) \)

where \( \tau \) is the tensor part of \(\chi\) that is the (inverse) metric by default

Parameters

lambda	input prefactor
phi	the vector to take the variation of
sigma	(out) the variation

Template Parameters

ContainerTypes

must be usable with Container in The dg dispatch system

variation() [3/3]

template<class Geometry , class Matrix , class Container >

template<class ContainerTypeL , class ContainerType0 , class ContainerType1 >

void dg::Elliptic3d< Geometry, Matrix, Container >::variation ( value_type  alpha, const ContainerTypeL &  lambda, const ContainerType0 &  phi, value_type  beta, ContainerType1 &  sigma  )

inline

\( \sigma = \alpha \lambda^2 (\nabla\phi\cdot\tau\cdot\nabla \phi) + \beta \sigma\)

where \( \tau \) is the tensor part of \(\chi\) that is the (inverse) metric by default

Parameters

alpha	scalar input prefactor
lambda	input prefactor
phi	the vector to take the variation of
beta	the output prefactor
sigma	(inout) the variation

Template Parameters

ContainerTypes

must be usable with Container in The dg dispatch system

weights()

template<class Geometry , class Matrix , class Container >

const Container & dg::Elliptic3d< Geometry, Matrix, Container >::weights ( ) const

inline

Return the weights making the operator self-adjoint.

i.e. the volume form

Returns

volume form including weights

---

The documentation for this class was generated from the following file:

• elliptic.h