How partial differential equations are specified in Montjoie

We describe here how equations are specified in Montjoie for 2-D and 3-D linear equations. We distinguish equations for which we can use H1 approximation (continuous finite elements), L2 approximation (Discontinuous Galerkin) or H(curl) approximation (so-called edge elements). The choice between continuous or discontinous approximation has to be considered carefully. While discontinuous Galerkin formulations allow to obtain more reliable solutions with a control of spurious modes, these formulations induce less efficient computations (in terms of computational time and storage). H(curl) approximation is currently only used for Maxwell's equations (and variants). When you implement a new set of equations, it is important to use it with compatible finite elements (for example, if you use a Discontinuous Galerkin formulation, a discontinuous finite element has to be used).

In Montjoie, we consider that the main unknown U (which can contain several components) is defined in all the computational domain and is solution of the same partial differential equation in this domain. Therefore, it is not "natural" to consider coupled problems with different partial differential equations in each part of the computational domain, and for which we could have considered different approximations (for example H(curl) approximation on a part of the domain and H1 in another part). In order to perform such operations, it is needed to introduce two objects EllipticProblem (or HyperbolicProblem for evolution problems), each object being devoted to the resolution of a different set of partial differential equations. Such things are performed for fluid-structure problems and vibroacoustic problems, but there is no general class developed for coupled problems.

General formulation used for H1 approximation

We are considering the following partial differential equations :

The unknown U can be scalar (i.e containing only one component) or vectorial. The tensors A(x), C(x), D(x) and E(x) can be constant or depending on space x, and on time t.

Variational formulation and how to specify a new set of equations

The variational formulation for continuous finite elements is the following one

where operator B(x) is describing the boundary condition imposed on Γ (Dirichlet-to-Neumann operator). We can also consider more complicated boundary conditions (cf. the section devoted to Boundary conditions). If we have D=E* and if A and C are symmetric, the variational formulation is symmetric even if the numerical integration is not exact. If you to specify tensors A, C, D and E, it is needed to construct a class with the following prototype :

//! inheritance from GenericEquation
//! 1 : scalar elements, 2 : edge elements, 3 : facet elements
class MyNewEquation
 : public GenericEquation<Complex_wp, 1 > 
{
  public :
    enum{nb_unknowns = 2, nb_components_en = 1, nb_components_hn = 1, nb_unknowns_scal = 2, nb_unknowns_vec = 4, nb_unknowns_hdg = 0};
    
    static const bool FormulationDG = false;
    static const bool FirstOrderFormulation = false;
    static const bool DiscontinuousDiMatrix = false;
    static const bool TensorStiffnessSymmetric = true;
  
    // computation of dV = C*dU
    template<class TypeElement, class TypeEquation, class Vector1, class Vector2>
    static void ApplyTensorStiffness(const EllipticProblem<TypeElement,TypeEquation>& var, int i, int j,
				     const NatureMatrix& nat_mat, int ref, Vector1& dU, Vector2& dV);

    // computation of V = A*U
    template<class TypeElement, class TypeEquation, class Vector1, class Vector2>
    static void ApplyTensorMass(const EllipticProblem<TypeElement,TypeEquation>& var, int i, int j,
				const NatureMatrix& nat_mat, int ref, Vector1& U, Vector2& V);

    // computation of V = E*dU
    template<class TypeElement, class TypeEquation, class Vector1, class Vector2>
    static void ApplyGradientUnknown(const EllipticProblem<TypeElement,TypeEquation>& var, int i, int j,
				     const NatureMatrix& nat_mat, int ref, Vector1& dU, Vector2& V);

    // computation of dV = D*U
    template<class TypeElement, class TypeEquation, class Vector1, class Vector2>
    static void ApplyGradientFctTest(const EllipticProblem<TypeElement,TypeEquation>& var, int i, int j,
				     const NatureMatrix& nat_mat, int ref, Vector1& U, Vector2& dV);

    // providing C in Cgrad_grad
    template<class TypeElement, class TypeEquation, class MatStiff>
    static void GetGradGradTensor(const EllipticProblem<TypeElement,TypeEquation>& vars,
                                  int i, int j, const NatureMatrix& nat_mat, int ref, MatStiff& Cgrad_grad);

    // providing A in mass
    template<class TypeElement, class TypeEquation, class Matrix1>
    static void GetTensorMass(const EllipticProblem<TypeElement, TypeEquation>& vars,
			      int i, int j, const NatureMatrix& nat_mat, int ref, Matrix1& mass);

    // providing D and E in Dgrad_phi and Ephi_grad
    template<class TypeElement, class TypeEquation, class MatStiff>
    static void GetGradPhiTensor(const EllipticProblem<TypeElement,TypeEquation>& vars, int i, int j,
                                 const NatureMatrix& nat_mat, int ref, MatStiff& Dgrad_phi, MatStiff& Ephi_grad);
};

As you can see, you need to specify both the tensors A, C, D, and E and also their product with "vectors". The methods requiring products are actually used for matrix-vector product (if the matrix is not stored) whereas the methods requiring the tensors are used to compute the elementary matrices. Since these functions are called a lot of times, it is important to optimize them by using precomputed arrays if needed, inlining them if they are placed in cxx files, and unrolling loops. FormulationDG should be set to false if you want to use continuous finite elements, true for discontinuous Galerkin. FirstOrderFormulation should be set to true if tensor C(x) is null (in which case you have only first-order derivatives in space). DiscontinuousDiMatrix and TensorStiffnessSymmetric are parameters used only for parallel execution, since continuous tensor D(x) and symmetric tensor C(x) are inducing lighter algorithms. For constant parameters (which are set once in the class, and are known at compilation) defined in the enum, we have the following convention :

The two last numbers are very specific to evolution problems, which are put into a first order problem through a mixed formulation. The constants nb_unknowns_en, nb_unknowns_hn are used when computing U and dU/dn on surfacic meshes (in order to write the correct number of components). The only important attribute is nb_unknowns. Details and explanations of methods in class GenericEquation are available below :

ApplyTensorStiffness Applies tensor C to a matrix dU
ApplyTensorMass Applies tensor A to a vector U
ApplyGradientUnknown Applies tensor E to a matrix dU
ApplyGradientFctTest Applies tensor D to a vector U
GetGradGradTensor fills tensor C
GetTensorMass fills tensor A
GetGradPhiTensor fills tensors D and E

One you have defined a class for your equation, you add a branch for this new problem with the following code :

// you can derive from VarHarmonic
// but if your problem is close to another one (like Helmholtz or Maxwell)
// you could derive from this closer problem
template<class TypeElement, class TypeEquation>
class MyNewProblem : public VarHarmonic<TypeElement, TypeEquation>
{
public :
// you can add attributes specific to the new equation
// usually the physical indices
VectReal_wp ref_lambda, ref_mu;

// you overload the methods you want
// for physical indices, usually we overload InitIndices and SetIndices
void InitIndices(int n);
void SetIndices(int i, const Vector<string>& param);

};

// then you create your specialized class EllipticProblem :
template<class TypeElement>
class EllipticProblem<TypeElement, MyNewEquation> : public MyNewProblem<TypeElement, TypeEquation>
{
  // nothing to put here usually
};

Then in the executable file (.cc), you can use this new equation :

int main(int argc, char** argv)
{
  InitMontjoie(argc, argv);

  if (argc != 2)
    {
      cout << "Give a data file " << endl;
      abort();
    }

  string input_file(argv[1]);
  EllipticProblem<TriangleLobatto, MyNewEquation> vars;
  vars.RunAll(input_file);

}

Time-harmonic/Time-domain computations

The tensors A, C, D and E can depend on the space (variable x), but they may also involve time derivatives, for instance :

In time-harmonic domain, it would be transformed into :

For time-harmonic computations, the tensors A(x), C(x), D(x) and E(x) are directly filled by using the pulsation (methods GetOmega(), GetSquareOmega() and GetMiomega()). For time-domain computations, we are only considering first and second-order derivatives in time. When the equations are only involving first derivatives (i.e. M(x) = 0), we say that the problem is of first-order. For implicit schemes, tensors C, D and E can contain derivatives in time, whereas explicit schemes are valid if only tensor A is involving derivatives in time. For implicit schemes, we need an intermediary object for the resolution of a stationary-like problem when tensor A would be replaced by :

Coefficients α, β and γ are given by the time-scheme used. For example, for a Crank-Nicolson scheme (M(x) = 0), we would have

The linear system to solve involves only Un+1, so we would have to solve :

These coefficients α, β and γ are contained in object "Nature_Matrix", through methods SetCoefMass/GetCoefMass (for α), SetCoefDamping/GetCoefDamping (for β) and SetCoefStiffness/GetCoefStiffness (for γ). For example, we show in the example below how to compute a matrix with different values of these coefficients.

// considering you have constructed an EllipticProblem object
void MyFunction(const EllipticProblem<TypeElement, TypeEquation>& var)
{
  // then you declare a "Nature_Matrix" object
  typename TypeEquation::Nature_Matrix nat_mat;
  
  // you set coefficients contained in this object :
  Real_wp dt = 0.01;
  nat_mat.SetCoefMass(0);
  nat_mat.SetCoefDamping(1.0/dt);
  nat_mat.SetCoefStiffness(0.5);

  // then you compute the finite element matrix with these coefficients
  Matrix<Real_wp, General, ArrayRowSparse> A;
  var.AddMatrixWithBC(A, nat_mat);

  // matrix is computed, you can factorize it if you want
  All_MatrixLU<Real_wp> mat_lu;
  mat_lu.Factorize(A);

}

As a result, the specification of time-domain equations is containing the type of the associated stationary equation. There are also two other methods GetMassCoefficient/GetDampingCoefficient, which are used in discontinuous Galerkin methods for explicit schemes. Here an example of a new time-domain equation

 
  class MyNewTimeEquation
  {
  public :
    // associated stationary equation
    typedef MyNewEquation TypeEquationStationary
    // using real numbers
    typedef Real_wp Complexe;
    
    //! mass coefficient
    template<class TypeElement, class TypeEquation>
    static inline Real_wp GetMassCoefficient(const HyperbolicProblem<TypeElement,TypeEquation>& var,
					     int i, int j, int m, int ref)
    { 
      return Real_wp(1);
    }
        
    //! damping coefficient
    template<class TypeElement, class TypeEquation>
    static inline Real_wp GetDampingCoefficient(const HyperbolicProblem<TypeElement,TypeEquation<& var,
						int i, int j, int m, int ref)
    { 
      return Real_wp(0);
    }        
  };

Unsteady problems are then derived with HyperbolicProblem class :

// usually it is derived from VarInstationary, but you could also
// derive from a closer class 
template<class TypeElement, class TypeEquation>
class MyNewTimeProblem : public VarInstationary<TypeElement, TypeEquation>
{
  public :

  // you may overload methods specific for this problem
  // usually nothing to do if only the equation has changed

};

// and finally the HyperbolicProblem class :
template<class TypeElement>
class HyperbolicProblem<TypeElement, MyNewTimeEquation>
 : public MyNewTimeProblem<TypeElement, MyNewTimeEquation>
{
  // nothing to do
};

In the main file .cc, the object HyperbolicProblem is then used to launch the simulation.

General formulation used for Discontinuous Galerkin approximation

We are considering the same partial differential equations as for continuous approximation :

The variational formulation is classical :

where ν is the normale outward from K- to K+, and we denote

The tensor denoted P is a penalization matrix. For first-order problems (C = 0), if you set :

Then the obtained LDG formulation is the classical one with upwind fluxes. You can also choose simpler forms of P, for example :

The specification of matrix P is done in methods GetPenalDG/MltPenalDG. When second-order formulation is chosen, matrix N will be equal to 0, and the boundary condition will be set will the same method as for continuous finite elements (through a DtN operator). When first-order formulation is set, the matrix N will contain the boundary condition operator. On each face of the boundary, u- is known and used in fluxes, while for u+, we are replacing the term :

As a result, you have to exploit boundary conditions to express Nb u+ as a function of u-. Matrix N is given through methods MltNabc/GetNabc. Finally, we show the header of a class implementing the resolution of a first-order problem with Discontinuous Galerkin method:

//! inheritance from GenericEquation
class MyNewEquation
 : public GenericEquation<Complex_wp, Symmetric, General, GlobalMatrixH1<Complex_wp> > 
{
  public :
    enum{nb_unknowns = 2, nb_components_en = 1, nb_components_hn = 1, nb_unknowns_scal = 2, nb_unknowns_vec = 4};
    
    static const bool FormulationDG = true;
    static const bool FirstOrderFormulation = true;
    static const bool DiscontinuousDiMatrix = false;
    static const bool TensorStiffnessSymmetric = true;
  
    // computation of dV = C*dU
    template
    static void ApplyTensorStiffness(const EllipticProblem<TypeElement,TypeEquation>& var, int i, int j,
				     const NatureMatrix& nat_mat, int ref, Vector1& dU, Vector2& dV);

    // computation of V = A*U
    template
    static void ApplyTensorMass(const EllipticProblem<TypeElement,TypeEquation>& var, int i, int j,
				const NatureMatrix& nat_mat, int ref, Vector1& U, Vector2& V);

    // computation of V = E*dU
    template
    static void ApplyGradientUnknown(const EllipticProblem<TypeElement,TypeEquation>& var, int i, int j,
				     const NatureMatrix& nat_mat, int ref, Vector1& dU, Vector2& V);

    // computation of dV = D*U
    template
    static void ApplyGradientFctTest(const EllipticProblem<TypeElement,TypeEquation>& var, int i, int j,
				     const NatureMatrix& nat_mat, int ref, Vector1& U, Vector2& dV);

    // providing C in Cgrad_grad
    template
    static void GetGradGradTensor(const EllipticProblem<TypeElement,TypeEquation>& vars,
                                  int i, int j, const NatureMatrix& nat_mat, int ref, MatStiff& Cgrad_grad);

    // providing A in mass
    template
    static void GetTensorMass(const EllipticProblem<TypeElement, TypeEquation>& vars,
			      int i, int j, const NatureMatrix& nat_mat, int ref, Matrix1& mass);

    // providing D and E in Dgrad_phi and Ephi_grad
    template
    static void GetGradPhiTensor(const EllipticProblem<TypeElement,TypeEquation>& vars, int i, int j,
                                 const NatureMatrix& nat_mat, int ref, MatStiff& Dgrad_phi, MatStiff& Ephi_grad);

    // providing N in Nabc
    template
    static void GetNabc(Matrix1& Nabc, typename Dimension::R_N& normale,
			int ref, int iquad, int k,const NatureMatrix& nat_mat, int ref_d,
			const GenericPb& vars, const TypeElt& Fb);

    // computation of Un = N*Vn
    template
    static void MltNabc(typename TypeElement::Dimension::R_N& normale, int ref,
                        const Vector1& Vn, Vector1& Un, int num_elem1, int k,
                        const NatureMatrix& nat_mat, int ref_d,
                        const EllipticProblem<TypeElement,TypeEquation>& vars, const TypeElt& Fb);
    
    // providing P in Nabc
    template
    static void GetPenalDG(Matrix1& Nabc, typename Dimension::R_N& normale,
			    int iquad, int k, const NatureMatrix& nat_mat, int ref, 
			    const GenericPb& vars, const TypeElt& Fb);
    
    // computation of Un = P*Vn
    template
    static void MltPenalDG(const typename Dimension::R_N& normale, const Vector1& Vn, Vector2& Un,
			    int i, int k, const NatureMatrix& nat_mat, int ref,
			    const GenericPb& vars, const TypeElt& Fb);
};

General formulation used for edge elements

Edge elements are used to approximate H(curl) finite element space. In Montjoie, they are only applied to Maxwell's equations (and variants). As a result, the variational formulation considered is quite simple

The associated pde is equal to :

Boundary condition is set with a DtN operator (denoted B(x) here). If you want to see how to implement tensors A(x) and C(x), you can look at classes HarmonicMaxwellEquation_2D and HarmonicMaxwellEquation_3D.