Elliptic operators

Elliptic \( -\nabla\cdot (\chi \nabla f)\) and Helmholtz \( (\chi + \alpha \Delta) f\). More...

Collaboration diagram for Elliptic operators:

Classes
class	dg::Elliptic1d< Geometry, Matrix, Container >
	A 1d negative elliptic differential operator \( -\partial_x ( \chi \partial_x ) \). More...
class	dg::Elliptic2d< Geometry, Matrix, Container >
	A 2d negative elliptic differential operator \( -\nabla \cdot ( \mathbf{\chi}\cdot \nabla ) \). More...
class	dg::Elliptic3d< Geometry, Matrix, Container >
	A 3d negative elliptic differential operator \( -\nabla \cdot ( \mathbf{\chi}\cdot \nabla ) \). More...
struct	dg::GeneralHelmholtz< Matrix, Container >
	A general Helmholtz-type operator \( (\chi-\alpha F) \). More...
struct	dg::Helmholtz2< Geometry, Matrix, Container >
	DEPRECATED, Matrix class that represents a more general Helmholtz-type operator. More...
class	dg::RefinedElliptic< Geometry, IMatrix, Matrix, Container >
	The refined version of Elliptic. More...

Typedefs
template<class Geometry , class Matrix , class Container >
using	dg::Elliptic = Elliptic2d<Geometry, Matrix, Container>
	A 2d negative elliptic differential operator \( -\nabla \cdot ( \mathbf{\chi}\cdot \nabla ) \).
template<class Geometry , class Matrix , class Container >
using	dg::Helmholtz = GeneralHelmholtz<dg::Elliptic2d<Geometry,Matrix,Container>, Container>
	a 2d Helmholtz opereator \( (\chi - \alpha F)\) with \( F = -\Delta\)
template<class Geometry , class Matrix , class Container >
using	dg::Helmholtz1d = GeneralHelmholtz<dg::Elliptic1d<Geometry,Matrix,Container>, Container>
	a 1d Helmholtz opereator \( (\chi - \alpha F)\) with \( F = -\partial_x^2\)
template<class Geometry , class Matrix , class Container >
using	dg::Helmholtz2d = GeneralHelmholtz<dg::Elliptic2d<Geometry,Matrix,Container>, Container>
	a 2d Helmholtz opereator \( (\chi - \alpha F)\) with \( F = -\Delta\)
template<class Geometry , class Matrix , class Container >
using	dg::Helmholtz3d = GeneralHelmholtz<dg::Elliptic3d<Geometry,Matrix,Container>, Container>
	a 3d Helmholtz opereator \( (\chi - \alpha F)\) with \( F = -\Delta\)

Detailed Description

Elliptic \( -\nabla\cdot (\chi \nabla f)\) and Helmholtz \( (\chi + \alpha \Delta) f\).

Typedef Documentation

Elliptic

template<class Geometry , class Matrix , class Container >

using dg::Elliptic = Elliptic2d<Geometry, Matrix, Container>

A 2d 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 \right)\right) + \partial_y\left(\sqrt{g} \left(\chi^{yx}\partial_x + \chi^{yy}\partial_y \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 Elliptic in an inversion problem

    //create an Elliptic object
    dg::Elliptic<dg::x::CartesianGrid2d, dg::x::DMatrix, dg::x::DVec>
        pol_forward( grid, dg::forward, jfactor);
 
    //Set the chi function (chi is a dg::x::DVec of size grid.size())
    pol_forward.set_chi( chi);
 
    //construct an pcg object
    dg::PCG pcg( x, n*n*Nx*Ny);
 
    //invert the elliptic equation
    pcg.solve( pol_forward, x, b, chi_inv, w2d, eps);
 
    //compute the error (solution contains analytic solution
    dg::blas1::axpby( 1.,x,-1., solution, error);
 
    //compute the L2 norm of the error
    double err = dg::blas2::dot( w2d, error);
 
    //output the relative error
    DG_RANK0 std::cout << " "<<sqrt( err/norm) << "\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

Helmholtz

template<class Geometry , class Matrix , class Container >

using dg::Helmholtz = GeneralHelmholtz<dg::Elliptic2d<Geometry,Matrix,Container>, Container>

a 2d Helmholtz opereator \( (\chi - \alpha F)\) with \( F = -\Delta\)

where \( \chi\) is a vector and \(\alpha\) a scalar and \(F\) is an operator.

Attention

If \( F\) is the Elliptic operator then the Laplacian in this formula becomes positive as opposed to the negative sign in the Elliptic operator

Can be used by the dg::PCG class. The following example shows how the class can be used to act as a Helmholtz2 operator:

    std::cout << "Alternative test with two Helmholtz operators\n";
    dg::Helmholtz< dg::CartesianGrid2d, dg::DMatrix, dg::DVec > gamma1inv(alpha, {grid2d, dg::centered});
    gamma1inv.set_chi( chi);
    dg::PCG<dg::DVec> pcgO(  x, grid2d.size());
    dg::PCG<dg::DVec> pcgOO( x, grid2d.size());
    t.tic();
    unsigned number1 = pcgO.solve( gamma1inv, phi, rholap, 1., w2d, eps/100);
    dg::blas1::pointwiseDot( phi, chi, phi);
    unsigned number2 = pcgOO.solve( gamma1inv, x, phi, 1., w2d, eps/100);
    t.toc();
    //Evaluation
    dg::blas1::axpby( 1., sol, -1., x);

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

Note

The intention is for Matrix to be one of the Elliptic classes with the weights() and precond() methods defined. If Matrix is to be an arbitrary functor then it is more convenient to simply directly use a lambda function to achieve the computation of \( y = (\chi - \alpha F)x\)

Helmholtz1d

template<class Geometry , class Matrix , class Container >

using dg::Helmholtz1d = GeneralHelmholtz<dg::Elliptic1d<Geometry,Matrix,Container>, Container>

a 1d Helmholtz opereator \( (\chi - \alpha F)\) with \( F = -\partial_x^2\)

where \( \chi\) is a vector and \(\alpha\) a scalar and \(F\) is an operator.

Attention

If \( F\) is the Elliptic operator then the Laplacian in this formula becomes positive as opposed to the negative sign in the Elliptic operator

Can be used by the dg::PCG class. The following example shows how the class can be used to act as a Helmholtz2 operator:

    std::cout << "Alternative test with two Helmholtz operators\n";
    dg::Helmholtz< dg::CartesianGrid2d, dg::DMatrix, dg::DVec > gamma1inv(alpha, {grid2d, dg::centered});
    gamma1inv.set_chi( chi);
    dg::PCG<dg::DVec> pcgO(  x, grid2d.size());
    dg::PCG<dg::DVec> pcgOO( x, grid2d.size());
    t.tic();
    unsigned number1 = pcgO.solve( gamma1inv, phi, rholap, 1., w2d, eps/100);
    dg::blas1::pointwiseDot( phi, chi, phi);
    unsigned number2 = pcgOO.solve( gamma1inv, x, phi, 1., w2d, eps/100);
    t.toc();
    //Evaluation
    dg::blas1::axpby( 1., sol, -1., x);

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

Note

The intention is for Matrix to be one of the Elliptic classes with the weights() and precond() methods defined. If Matrix is to be an arbitrary functor then it is more convenient to simply directly use a lambda function to achieve the computation of \( y = (\chi - \alpha F)x\)

Helmholtz2d

template<class Geometry , class Matrix , class Container >

using dg::Helmholtz2d = GeneralHelmholtz<dg::Elliptic2d<Geometry,Matrix,Container>, Container>

a 2d Helmholtz opereator \( (\chi - \alpha F)\) with \( F = -\Delta\)

where \( \chi\) is a vector and \(\alpha\) a scalar and \(F\) is an operator.

Attention

If \( F\) is the Elliptic operator then the Laplacian in this formula becomes positive as opposed to the negative sign in the Elliptic operator

Can be used by the dg::PCG class. The following example shows how the class can be used to act as a Helmholtz2 operator:

    std::cout << "Alternative test with two Helmholtz operators\n";
    dg::Helmholtz< dg::CartesianGrid2d, dg::DMatrix, dg::DVec > gamma1inv(alpha, {grid2d, dg::centered});
    gamma1inv.set_chi( chi);
    dg::PCG<dg::DVec> pcgO(  x, grid2d.size());
    dg::PCG<dg::DVec> pcgOO( x, grid2d.size());
    t.tic();
    unsigned number1 = pcgO.solve( gamma1inv, phi, rholap, 1., w2d, eps/100);
    dg::blas1::pointwiseDot( phi, chi, phi);
    unsigned number2 = pcgOO.solve( gamma1inv, x, phi, 1., w2d, eps/100);
    t.toc();
    //Evaluation
    dg::blas1::axpby( 1., sol, -1., x);

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

Note

The intention is for Matrix to be one of the Elliptic classes with the weights() and precond() methods defined. If Matrix is to be an arbitrary functor then it is more convenient to simply directly use a lambda function to achieve the computation of \( y = (\chi - \alpha F)x\)

Helmholtz3d

template<class Geometry , class Matrix , class Container >

using dg::Helmholtz3d = GeneralHelmholtz<dg::Elliptic3d<Geometry,Matrix,Container>, Container>

a 3d Helmholtz opereator \( (\chi - \alpha F)\) with \( F = -\Delta\)

where \( \chi\) is a vector and \(\alpha\) a scalar and \(F\) is an operator.

Attention

If \( F\) is the Elliptic operator then the Laplacian in this formula becomes positive as opposed to the negative sign in the Elliptic operator

Can be used by the dg::PCG class. The following example shows how the class can be used to act as a Helmholtz2 operator:

    std::cout << "Alternative test with two Helmholtz operators\n";
    dg::Helmholtz< dg::CartesianGrid2d, dg::DMatrix, dg::DVec > gamma1inv(alpha, {grid2d, dg::centered});
    gamma1inv.set_chi( chi);
    dg::PCG<dg::DVec> pcgO(  x, grid2d.size());
    dg::PCG<dg::DVec> pcgOO( x, grid2d.size());
    t.tic();
    unsigned number1 = pcgO.solve( gamma1inv, phi, rholap, 1., w2d, eps/100);
    dg::blas1::pointwiseDot( phi, chi, phi);
    unsigned number2 = pcgOO.solve( gamma1inv, x, phi, 1., w2d, eps/100);
    t.toc();
    //Evaluation
    dg::blas1::axpby( 1., sol, -1., x);

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

Note

The intention is for Matrix to be one of the Elliptic classes with the weights() and precond() methods defined. If Matrix is to be an arbitrary functor then it is more convenient to simply directly use a lambda function to achieve the computation of \( y = (\chi - \alpha F)x\)