Discontinuous Galerkin Library
#include "dg/algorithm.h"
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 x = dg::evaluate( initial, grid);
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";
A 3d negative elliptic differential operator .
Definition: elliptic.h:545
Preconditioned conjugate gradient method to solve .
Definition: pcg.h:57
@ centered
centered derivative (cell to the left and right and current cell)
Definition: enums.h:100
@ x
x direction
thrust::host_vector< real_type > evaluate(UnaryOp f, const RealGrid1d< real_type > &g)
Evaluate a 1d function on grid coordinates.
Definition: evaluation.h:67
get_host_vector< Geometry > volume(const Geometry &g)
Create the volume element on the grid (including weights!!)
Definition: transform.h:225
thrust::device_vector< double > DVec
Device Vector. The device can be an OpenMP parallelized cpu or a gpu. This depends on the value of th...
Definition: typedefs.h:23
three-dimensional Grid with Cylindrical metric
Definition: base_geometry.h:271
Template Parameters
GeometryA 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:
MatrixA 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:
ContainerA 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
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
gThe Grid; boundary conditions are taken from here
dirDirection 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_jumpIf 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
gThe Grid
bcxboundary condition in x
bcyboundary contition in y
bczboundary contition in z
dirDirection 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_jumpIf 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
Paramsdeduced by the compiler
Parameters
psparameters 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
sigmaThe 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
ContainerType0must 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
tauThe new tensor part in \(\chi\) (must be positive definite)
Note
the 3d parts in tau will be ignored for 2d computations
Template Parameters
ContainerType0must 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_2dif 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_jfactorThe 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_weightingSwitch 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
xleft-hand-side
yresult
Template Parameters
ContainerTypesmust 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
alphaa scalar
xleft-hand-side
betaa scalar
yresult
Template Parameters
ContainerTypesmust 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
phithe vector to take the variation of
sigma(inout) the variation
Template Parameters
ContainerTypesmust 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
lambdainput prefactor
phithe vector to take the variation of
sigma(out) the variation
Template Parameters
ContainerTypesmust 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
alphascalar input prefactor
lambdainput prefactor
phithe vector to take the variation of
betathe output prefactor
sigma(inout) the variation
Template Parameters
ContainerTypesmust 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: