Solving Partial Differential Equations

1. Introduction

This tutorial demonstrates how to solve partial differential equations (PDEs) using Feel++'s finite element framework. We’ll cover the complete workflow from problem formulation to solution visualization.

2. The Finite Element Method Workflow

Solving a PDE with Feel++ follows these steps:

Problem Formulation: Define the weak form of your PDE
Mesh Generation: Create or load a computational domain
Function Space Creation: Define trial and test spaces
Form Assembly: Build bilinear and linear forms
Boundary Conditions: Apply essential and natural BCs
System Solution: Solve the resulting linear/nonlinear system
Post-processing: Analyze and visualize results

3. Example 1: Poisson Equation

Let’s solve the Poisson equation as our first example:

-Δu = f    in Ω
u = g      on ∂Ω

3.1. Weak Formulation

Find u ∈ V such that: ∫_Ω ∇u·∇v dx = ∫_Ω fv dx ∀v ∈ V₀

3.2. Implementation

#include <feel/feel.hpp>

int main( int argc, char** argv )
{
    using namespace Feel;
    using namespace Feel::vf;

    // Define command line options
    po::options_description laplacian_options( "Laplacian solver options" );
    laplacian_options.add_options()
        ( "order", po::value<int>()->default_value( 2 ), "polynomial order" )
        ( "penalbc", po::value<double>()->default_value( 100 ), "penalization parameter for Dirichlet BC" )
        ( "weak-bc", po::value<bool>()->default_value( false ), "use weak boundary conditions" )
        ( "mu", po::value<double>()->default_value( 1.0 ), "diffusion coefficient" );

    Environment env( _argc=argc, _argv=argv,
                     _desc=laplacian_options,
                     _about=about( _name="laplacian",
                                   _author="Feel++ Consortium",
                                   _email="feelpp-devel@feelpp.org" ) );

    // Load mesh
    auto mesh = loadMesh(_mesh=new Mesh<Simplex<2>>);
    
    // Create function space
    int order = ioption("order");
    auto Vh = Pch<2>( mesh );
    auto u = Vh->element("u");
    auto v = Vh->element("v");
    
    // Get parameters
    double mu = doption("mu");
    double penalbc = doption("penalbc");
    bool weak_bc = boption("weak-bc");
    
    LOG(INFO) << "Solving Laplacian equation with:";
    LOG(INFO) << "  - polynomial order: " << order;
    LOG(INFO) << "  - diffusion coefficient: " << mu;
    LOG(INFO) << "  - weak BC: " << (weak_bc ? "yes" : "no");
    
    // Define exact solution and right-hand side
    auto u_exact = sin(pi*Px())*sin(pi*Py());
    auto f = mu * 2*pi*pi*sin(pi*Px())*sin(pi*Py());  // -div(mu*grad(u)) = f
    auto g = u_exact;  // Dirichlet boundary condition
    
    // Assemble the bilinear form a(u,v) = ∫ μ ∇u·∇v
    auto a = form2( _trial=Vh, _test=Vh );
    a = integrate( _range=elements(mesh),
                   _expr=mu * gradt(u) * trans(grad(v)) );
    
    // Add boundary conditions
    if ( weak_bc )
    {
        // Weak Dirichlet boundary conditions using Nitsche's method
        a += integrate( _range=boundaryfaces(mesh),
                       _expr=-mu * gradt(u) * N() * id(v)
                            -mu * grad(v) * N() * idt(u)
                            + penalbc * idt(u) * id(v) / hFace() );
    }
    else
    {
        // Strong Dirichlet boundary conditions
        a += on( _range=boundaryfaces(mesh), _rhs=l, _element=u, _expr=g );
    }
    
    // Assemble the linear form l(v) = ∫ f·v
    auto l = form1( _test=Vh );
    l = integrate( _range=elements(mesh),
                   _expr=f * id(v) );
    
    if ( weak_bc )
    {
        // Add boundary terms for weak BC
        l += integrate( _range=boundaryfaces(mesh),
                       _expr=-mu * g * grad(v) * N() + penalbc * g * id(v) / hFace() );
    }
    
    // Solve the linear system
    a.solve( _rhs=l, _solution=u );
    
    // Compute error
    auto error_l2 = normL2( _range=elements(mesh), _expr=idv(u) - u_exact );
    auto error_h1 = normH1( _range=elements(mesh), 
                           _expr=idv(u) - u_exact,
                           _grad_expr=gradv(u) - grad<2>(u_exact) );
    
    LOG(INFO) << "L2 error: " << error_l2;
    LOG(INFO) << "H1 error: " << error_h1;
    
    // Export results
    auto e = exporter( _mesh=mesh, _name="laplacian" );
    e->add( "solution", u );
    e->add( "exact", vf::project(_space=Vh, _range=elements(mesh), _expr=u_exact) );
    e->add( "error", vf::project(_space=Vh, _range=elements(mesh), _expr=abs(idv(u) - u_exact)) );
    e->save();
    
    return 0;
}

3.3. Key Components Explained

Bilinear Form: a(u,v) = ∫ μ ∇u·∇v

auto a = form2( _trial=Vh, _test=Vh );
a = integrate( _range=elements(mesh),
               _expr=mu * gradt(u) * trans(grad(v)) );

Linear Form: l(v) = ∫ fv

auto l = form1( _test=Vh );
l = integrate( _range=elements(mesh),
               _expr=f * id(v) );

Boundary Conditions: Essential (Dirichlet) BC

a += on( _range=boundaryfaces(mesh), _rhs=l, _element=u, _expr=g );

4. Example 2: Stokes Equations

The Stokes equations model viscous flow:

-ν Δu + ∇p = f    in Ω  (momentum)
∇·u = 0           in Ω  (continuity)
u = g             on ∂Ω (no-slip)

4.1. Mixed Formulation

Find (u,p) ∈ V × Q such that: ∫_Ω ν ∇u:∇v dx + ∫_Ω p ∇·v dx = ∫_Ω f·v dx ∀v ∈ V₀ ∫_Ω q ∇·u dx = 0 ∀q ∈ Q

4.2. Implementation

#include <feel/feel.hpp>

int main( int argc, char** argv )
{
    using namespace Feel;
    using namespace Feel::vf;

    // Define options specific to Stokes problem
    po::options_description stokes_options( "Stokes solver options" );
    stokes_options.add_options()
        ( "nu", po::value<double>()->default_value( 1.0 ), "kinematic viscosity" )
        ( "penalbc", po::value<double>()->default_value( 100 ), "penalization for boundary conditions" );

    Environment env( _argc=argc, _argv=argv,
                     _desc=stokes_options,
                     _about=about( _name="stokes",
                                   _author="Feel++ Consortium",
                                   _email="feelpp-devel@feelpp.org" ) );

    // Load mesh
    auto mesh = loadMesh(_mesh=new Mesh<Simplex<2>>);
    
    // Create mixed function space: Velocity (P2) + Pressure (P1)
    auto Vh = Pchv<2>( mesh );     // Velocity space (vector-valued P2)
    auto Qh = Pch<1>( mesh );      // Pressure space (scalar P1) 
    auto Wh = Vh * Qh;             // Mixed space

    // Create elements
    auto U = Wh->element("(u,p)");
    auto u = U.template element<0>("u");  // velocity
    auto p = U.template element<1>("p");  // pressure
    auto V = Wh->element("(v,q)");
    auto v = V.template element<0>("v");  // test velocity
    auto q = V.template element<1>("q");  // test pressure
    
    // Problem parameters
    double nu = doption("nu");
    double penalbc = doption("penalbc");
    
    LOG(INFO) << "Solving Stokes equation with:";
    LOG(INFO) << "  - kinematic viscosity: " << nu;
    LOG(INFO) << "  - velocity DoFs: " << Vh->nDof();
    LOG(INFO) << "  - pressure DoFs: " << Qh->nDof();
    LOG(INFO) << "  - total DoFs: " << Wh->nDof();

    // Define manufactured solution for validation
    auto u1_exact = sin(pi*Px())*cos(pi*Py());
    auto u2_exact = -cos(pi*Px())*sin(pi*Py());
    auto p_exact = sin(pi*Px())*sin(pi*Py());
    
    // Compute corresponding body force
    auto f1 = -nu * (-2*pi*pi*sin(pi*Px())*cos(pi*Py())) + pi*cos(pi*Px())*sin(pi*Py());
    auto f2 = -nu * (2*pi*pi*cos(pi*Px())*sin(pi*Py())) + pi*sin(pi*Px())*cos(pi*Py());
    
    // Assemble Stokes system
    auto a = form2( _trial=Wh, _test=Wh );
    
    // Viscous term: ∫ ν ∇u : ∇v
    a += integrate( _range=elements(mesh),
                   _expr=nu * inner(gradt(u), grad(v)) );
    
    // Divergence of velocity (constraint): ∫ div(u) * q
    a += integrate( _range=elements(mesh),
                   _expr=divt(u) * id(q) );
    
    // Pressure gradient term: ∫ p * div(v)  
    a += integrate( _range=elements(mesh),
                   _expr=idt(p) * div(v) );
    
    // Right-hand side
    auto l = form1( _test=Wh );
    l += integrate( _range=elements(mesh),
                   _expr=f1*id(v).comp(X) + f2*id(v).comp(Y) );
    
    // Apply boundary conditions (no-slip walls)
    auto u_exact_vec = vec(u1_exact, u2_exact);
    a += on( _range=boundaryfaces(mesh), _rhs=l, _element=u, _expr=u_exact_vec );
    
    // Solve the system
    LOG(INFO) << "Assembling and solving Stokes system...";
    a.solve( _rhs=l, _solution=U );
    
    // Compute errors
    auto error_u_l2 = normL2( _range=elements(mesh), 
                             _expr=idv(u) - u_exact_vec );
    auto error_p_l2 = normL2( _range=elements(mesh),
                             _expr=idv(p) - p_exact );
    
    auto error_u_h1 = normH1( _range=elements(mesh),
                             _expr=idv(u) - u_exact_vec,
                             _grad_expr=gradv(u) - grad<2>(u_exact_vec) );
    
    LOG(INFO) << "Velocity L2 error: " << error_u_l2;
    LOG(INFO) << "Velocity H1 error: " << error_u_h1;
    LOG(INFO) << "Pressure L2 error: " << error_p_l2;
    
    // Check divergence-free condition
    auto div_error = normL2( _range=elements(mesh), _expr=divv(u) );
    LOG(INFO) << "Divergence error: " << div_error;
    
    // Export results
    auto e = exporter( _mesh=mesh, _name="stokes" );
    e->add( "velocity", u );
    e->add( "pressure", p );
    e->add( "velocity_exact", vf::project(_space=Vh, _range=elements(mesh), _expr=u_exact_vec) );
    e->add( "pressure_exact", vf::project(_space=Qh, _range=elements(mesh), _expr=p_exact) );
    e->add( "velocity_magnitude", vf::project(_space=Qh, _range=elements(mesh), _expr=sqrt(trans(idv(u))*idv(u))) );
    e->save();
    
    return 0;
}

5. Boundary Condition Types

Feel++ supports various boundary condition types:

5.1. Essential (Dirichlet) Boundary Conditions

Specify function values on boundaries:

// Strong enforcement
a += on( _range=boundaryfaces(mesh), _rhs=l, _element=u, _expr=g );

// Weak enforcement (Nitsche method)
a += integrate( _range=boundaryfaces(mesh),
               _expr=-gradt(u)*N()*id(v) - grad(v)*N()*idt(u)
                    + penalbc*idt(u)*id(v)/hFace() );
l += integrate( _range=boundaryfaces(mesh),
               _expr=-g*grad(v)*N() + penalbc*g*id(v)/hFace() );

5.2. Natural (Neumann) Boundary Conditions

Specify flux or normal derivatives:

// Natural BC: ∂u/∂n = h on boundary
l += integrate( _range=boundaryfaces(mesh),
               _expr=h * id(v) );

5.3. Robin Boundary Conditions

Mixed boundary conditions:

// Robin BC: αu + β∂u/∂n = γ
a += integrate( _range=boundaryfaces(mesh),
               _expr=alpha * idt(u) * id(v) );
l += integrate( _range=boundaryfaces(mesh),
               _expr=gamma * id(v) );

6. Advanced Solution Techniques

6.1. Nonlinear Problems

For nonlinear PDEs, use Newton’s method:

auto a = form2( _trial=Vh, _test=Vh );
auto l = form1( _test=Vh );
auto u = Vh->element();

// Newton iteration
for( int iter = 0; iter < max_iter; ++iter )
{
    // Assemble Jacobian and residual
    a.zero();
    l.zero();

    // Add linearized terms
    a += integrate( _range=elements(mesh),
                   _expr=/* linearized bilinear form */ );
    l += integrate( _range=elements(mesh),
                   _expr=/* residual form */ );

    // Apply BC and solve
    a += on( _range=boundaryfaces(mesh), _rhs=l, _element=du, _expr=cst(0.) );
    auto du = a.solve( _rhs=l, _solution=du );

    // Update solution
    u += du;

    // Check convergence
    if( du.l2Norm() < tolerance ) break;
}

6.2. Time-Dependent Problems

For time evolution, use time stepping schemes:

// BDF2 time scheme for ∂u/∂t - Δu = f
auto dt = doption("dt");
auto bdf2 = bdf( _space=Vh, _name="u" );

for( auto time : bdf2 )
{
    // Update time-dependent data
    auto f = expr(soption("f"), symbols, {{"t", time}});

    // Mass matrix contribution
    auto a = form2( _trial=Vh, _test=Vh );
    a = integrate( _range=elements(mesh),
                   _expr=bdf2.polyDerivCoefficient(0) * idt(u) * id(v) );

    // Stiffness matrix
    a += integrate( _range=elements(mesh),
                   _expr=gradt(u) * trans(grad(v)) );

    // Right-hand side
    auto l = form1( _test=Vh );
    l = integrate( _range=elements(mesh),
                   _expr=f * id(v) + bdf2.polyDeriv() * id(v) );

    // Solve and update
    a.solve( _rhs=l, _solution=u );
    bdf2.next(u);
}

7. Error Analysis and Validation

7.1. Manufactured Solutions

Test your implementation with known solutions:

// Define exact solution
auto u_exact = sin(pi*Px())*sin(pi*Py());

// Compute corresponding source term
auto f = 2*pi*pi*sin(pi*Px())*sin(pi*Py());

// Solve problem with this source
// ... problem setup and solve ...

// Compute errors
double error_l2 = normL2( _range=elements(mesh),
                         _expr=idv(u_h) - u_exact );
double error_h1 = normH1( _range=elements(mesh),
                         _expr=idv(u_h) - u_exact,
                         _grad_expr=gradv(u_h) - grad<2>(u_exact) );

7.2. Convergence Studies

Study how errors decrease with mesh refinement:

std::vector<double> h_values, l2_errors, h1_errors;

for( double h = 0.1; h >= 0.01; h /= 2 )
{
    // Generate mesh with characteristic size h
    auto mesh = createMesh( h );

    // Solve problem
    auto u_h = solve_problem( mesh );

    // Compute errors
    double error_l2 = compute_l2_error( u_h );
    double error_h1 = compute_h1_error( u_h );

    h_values.push_back( h );
    l2_errors.push_back( error_l2 );
    h1_errors.push_back( error_h1 );
}

// Compute convergence rates
for( int i = 1; i < h_values.size(); ++i )
{
    double rate_l2 = log(l2_errors[i]/l2_errors[i-1]) / log(h_values[i]/h_values[i-1]);
    double rate_h1 = log(h1_errors[i]/h1_errors[i-1]) / log(h_values[i]/h1_values[i-1]);

    LOG(INFO) << "h=" << h_values[i] << ", L2 rate=" << rate_l2 << ", H1 rate=" << rate_h1;
}

8. Solver Configuration

Feel++ provides extensive linear solver options:

8.1. Direct Solvers

For small to medium problems:

// In config file
ksp-type=preonly
pc-type=lu
pc-factor-mat-solver-package=mumps

8.2. Iterative Solvers

For large problems:

// Conjugate Gradient with AMG preconditioner
ksp-type=cg
pc-type=gamg
ksp-rtol=1e-8
ksp-max_it=1000

8.3. Problem-Specific Solvers

For Stokes problems:

// Block preconditioner for Stokes
pc-type=fieldsplit
pc-fieldsplit-type=schur
pc-fieldsplit-schur-factorization-type=diag
fieldsplit-0-ksp-type=preonly
fieldsplit-0-pc-type=ml
fieldsplit-1-ksp-type=preonly
fieldsplit-1-pc-type=ml

9. Performance Optimization

9.1. Memory Usage

Choose appropriate polynomial orders
Use discontinuous Galerkin for convection-dominated problems
Consider static condensation for high-order methods

9.2. Parallel Scaling

Use appropriate mesh partitioning
Configure solvers for your problem size
Profile communication patterns

9.3. Assembly Optimization

Vectorize expressions where possible
Use compile-time polynomial orders
Consider assembly time vs solution time trade-offs

10. Debugging and Troubleshooting

10.1. Common Issues

Matrix Singularity: - Check boundary conditions - Verify problem well-posedness - Inspect null space

Poor Convergence: - Improve preconditioners - Check mesh quality - Verify weak formulation

Incorrect Results: - Validate with manufactured solutions - Check unit consistency - Verify boundary condition implementation

10.2. Diagnostic Tools

// Check matrix properties
LOG(INFO) << "Matrix rows: " << a.matrixPtr()->size1();
LOG(INFO) << "Matrix cols: " << a.matrixPtr()->size2();
LOG(INFO) << "Matrix nnz: " << a.matrixPtr()->nnz();

// Examine solution properties
LOG(INFO) << "Solution L2 norm: " << u.l2Norm();
LOG(INFO) << "Solution min/max: " << u.min() << "/" << u.max();

11. Best Practices

Start Simple: Begin with known solutions for validation
Verify Convergence: Always check spatial and temporal convergence
Use Appropriate Orders: Match polynomial order to solution regularity
Monitor Performance: Profile assembly and solution times
Document Parameters: Keep track of solver and discretization parameters
Validate Physics: Ensure your formulation matches the physical problem

This comprehensive guide provides the foundation for solving PDEs with Feel++. The key is understanding the weak formulation, proper discretization, and effective solution strategies.