Mixed boundary conditions

[zhyxue0]

I try to solve the N-S equations, the velocity is the periodic boundary and the pressure is the Dirichlet boundary, how do I define the Functionspace?
Do I need to define two sets of mehes? mesh_1 = PeriodicRectangleMesh and mesh_2 = RectangleMesh

[dham]

I've converted this to a discussion since it seems to be a question about how to use Firedrake rather than an issue with Firedrake.

I think I don't understand your question. A periodic "boundary" is not a boundary at all, so it doesn't have boundary conditions. Could you please post the mathematical statement of the variational problem (or the strong form PDE) that you are trying to solve?

[zhyxue0]

Thanks, I summed up the problem.

Motivation

I want to do the accuracy tests, consider the following the incompressible Navier–Stokes equations

$\boldsymbol{u}_t-\nu \Delta \boldsymbol{u}+\boldsymbol{u} \cdot \nabla \boldsymbol{u}+\nabla p = \boldsymbol{0} $
$\nabla \cdot \boldsymbol{u} = 0$

in the $\Omega = [0, 2\pi]\times[0,2\pi]$

As the accuracy test, take the exact solution as

$\boldsymbol{u}_{\text{exact}}(x, y, t)=\left(e^{-2 t} \sin(x)\cos(y), -e^{-2 t} \cos (x) \sin (y)\right)$

$p_{\text{exact}}(x,y,t) = \frac{1}{4}e^{-4t}\left(\cos(2x) + \cos(2y)\right)$

It is easy to check that $\boldsymbol{u}_{\text{exact}}$ is divergence free ,

$\left.\boldsymbol{{u}{\text{exact}}}\right|{\partial \Omega}=0 \$

and $\int_{\Omega} p_{\text{exact}} dx = 0$

Temporal discretization

Consider the pressure-correction scheme:

Step 1: find $\boldsymbol{\tilde{u}} ^{n+1}$ by solving

$\frac{\boldsymbol{\tilde{u}} ^{n+1}- \boldsymbol{u}^{n}}{\Delta t} + \boldsymbol{u}^{n} \cdot \nabla \boldsymbol{\tilde{u}}^{n+1} - \nu \Delta \boldsymbol{\tilde{u}}^{n+1} + \nabla p^{n} = 0$

Step 2: find $(\boldsymbol{u}^{n+1}, p^{n+1})$ by solving

$\frac{\boldsymbol{u}^{n+1} - \boldsymbol{\tilde{u}}^{n+1}}{\Delta t} + \nabla ( p^{n+1} - p^{n}) = 0 $
$\nabla \cdot \boldsymbol{u}^{n+1} = 0$

Spatial discretization

LDG method
with alternating numerical fluxes

Step 1:

$\boldsymbol{\tilde{u}} ^{n+1} + \Delta t (\boldsymbol{u}^{n} \cdot \nabla )\boldsymbol{\tilde{u}}^{n+1} - \Delta t\nu \nabla \cdot \boldsymbol{Q}^{n+1} + \Delta t \boldsymbol{r}^{n} = \boldsymbol{u}^{n} $
$\boldsymbol{Q}^{n+1} - \nabla \boldsymbol{\tilde{u}}^{n+1} = 0 $
$\boldsymbol{r}^{n} = \nabla p^{n}$

Step 2:

$\boldsymbol{r}^{n+1} - \nabla p^{n+1} = 0$

$\boldsymbol{u}^{n+1} + \Delta t \boldsymbol{r}^{n+1} = \boldsymbol{\tilde{u}}^{n+1} + \Delta t \boldsymbol{r}^{n}$

$\nabla \cdot \boldsymbol{u}^{n+1} = 0$

Boundary conditions:

We consider different boundary conditions:

case 1

Periodic boundary conditions for both velocity $\boldsymbol{u}$ and pressure $p$, so add a constraint of $\int_{\Omega} p dx =0$

The treatment is similar to the Poisson equation with pure Neumann boundary condition. There are two options

1. Similar to the tutorial [Singular operators in mixed spaces](https://firedrakeproject.org/solving-interface.html#id21)[[¶](https://firedrakeproject.org/solving-interface.html#singular-operators-in-mixed-spaces)](https://firedrakeproject.org/solving-interface.html#singular-operators-in-mixed-spaces), call nullspace, but it doesn't seem to work. In step 2, pressure $p^{n+1}$ blows up.

2. The zero integral mean condition for the pressure $p$ is imposed via a real Lagrange multiplier $\lambda$, like in the tutorial The (R) space,

   Step 2' find $(\boldsymbol{u}^{n+1}, p^{n+1}, \lambda)$ by solving

   $\frac{\boldsymbol{u}^{n+1} - \boldsymbol{\tilde{u}}^{n+1}}{\Delta t} + \nabla ( p^{n+1} - p^{n}) = 0$

   $\int_{\Omega} \left(\nabla \cdot (\boldsymbol{u}^{n+1} - \boldsymbol{\tilde{u}}^{n+1} ) + \lambda \right)dx = 0$

​ $\int_{\Omega} p^{n+1} dx = 0$

Question

1. Why is nullspace not working properly?
2. R space corresponds to a dense matrix, but I don't know how to set up an efficient solver.

Here is my code

from firedrake import *
import math
import numpy as np

def Solver_Lagrange_multiplier(Mixed_space_2, un, r_n, u_tilde):
    r_trial, u_trial, p_trial, lambda_trial = TrialFunctions(Mixed_space_2)
    theta_3, theta_4, phi_1, mu = TestFunctions(Mixed_space_2)

    p_trial_flux = conditional(un('+') > 0, p_trial('-'), p_trial('+'))
    u_trial_flux = conditional(un('+') > 0, u_trial('+'), u_trial('-'))
    u_tilde_flux = conditional(un('+') > 0, u_tilde('+'), u_tilde('-'))

    # Eq.(2.1)'
    LHS_2 = dot(r_trial, theta_3)*dx \
        + p_trial * div(theta_3)*dx \
        - p_trial_flux * jump(theta_3, n)*dS 


    # Eq.(2.2)'
    LHS_2 += dot(u_trial, theta_4)*dx \
        + dtc*( dot(r_trial, theta_4)*dx)

    RHS_2 = dot(u_tilde, theta_4)*dx \
        + dtc*( dot(r_n, theta_4)*dx)

    # Eq.(2.3)'
    LHS_2 += - dot(u_trial, grad(phi_1))*dx \
            + dot(u_trial_flux, jump(phi_1,n))*dS 

    # add the constrain for p_{n+1}
    LHS_2 += lambda_trial * phi_1*dx \
        + p_trial*mu*dx    

    RHS_2 += - dot(u_tilde, grad(phi_1))*dx \
            + dot(u_tilde_flux, jump(phi_1,n))*dS 

    Solution_Mixed_space_2 = Function(Mixed_space_2)
    prob_2 = LinearVariationalProblem(LHS_2, RHS_2, Solution_Mixed_space_2)
    solv_2 = LinearVariationalSolver(prob_2)



    solv_2.solve()

    r = Function(W_h)
    lambda_1 = Function(R)
    r.assign(Solution_Mixed_space_2.sub(0))
    u.assign(Solution_Mixed_space_2.sub(1))
    p.assign(Solution_Mixed_space_2.sub(2))
    lambda_1.assign(Solution_Mixed_space_2.sub(3))

    print("R_space, u_max:", np.abs( Solution_Mixed_space_2.sub(1).vector().array() ).max())
    print("R_space, p_max:", np.abs( Solution_Mixed_space_2.sub(2).vector().array() ).max())


def Solver_nullspace(Mixed_space_2, un, r_n, u_tilde, nullmodes):
    
    r_trial, u_trial, p_trial = TrialFunctions(Mixed_space_2)
    theta_3, theta_4, phi_1 = TestFunctions(Mixed_space_2)

    p_trial_flux = conditional(un('+') > 0, p_trial('-'), p_trial('+'))
    u_trial_flux = conditional(un('+') > 0, u_trial('+'), u_trial('-'))

    # Eq.(2.1)
    LHS_2 = dot(r_trial, theta_3)*dx \
        + p_trial * div(theta_3)*dx \
        - p_trial_flux * jump(theta_3, n)*dS 


    # Eq.(2.2)
    LHS_2 += dot(u_trial, theta_4)*dx \
        + dtc*( dot(r_trial, theta_4)*dx)

    RHS_2 = dot(u_tilde, theta_4)*dx \
        + dtc*( dot(r_n, theta_4)*dx)

    # Eq.(2.3)
    LHS_2 += - dot(u_trial, grad(phi_1))*dx \
            + dot(u_trial_flux, jump(phi_1,n))*dS 

    Solution_Mixed_space_2 = Function(Mixed_space_2)
    prob_3 = LinearVariationalProblem(LHS_2, RHS_2, Solution_Mixed_space_2)

    # nullspace solver
    solv_3 = LinearVariationalSolver(prob_3, nullspace=nullmodes)
    solv_3.solve()

    print("nullspace solver, u_max:", np.abs( Solution_Mixed_space_2.sub(1).vector().array() ).max())
    print("nullspace solver, p_max:", np.abs( Solution_Mixed_space_2.sub(2).vector().array() ).max())

    # No nullspace
    solv_4 = LinearVariationalSolver(prob_3)
    solv_4.solve()

    print("No nullspace, u_max:", np.abs( Solution_Mixed_space_2.sub(1).vector().array() ).max())
    print("No nullspace, p_max:", np.abs( Solution_Mixed_space_2.sub(2).vector().array() ).max())


# define the mesh paramters
degree = 1
Nx = 20; Ny = Nx
Lx = 2*math.pi; Ly = 2*math.pi
Hx = Lx/Nx

# define the equation parameters
nu = 1.0

# define the time step
t = 0.0
dt = 1E-6
dtc = Constant(dt)
tc = Constant(t)

# define the mesh
mesh = PeriodicRectangleMesh(Nx, Ny, Lx, Ly, quadrilateral=False, reorder=None, diagonal="left")

# Mesh entities
n = FacetNormal(mesh)
h = CellDiameter(mesh)
x, y = SpatialCoordinate(mesh)

# define the orientation
Orientation = Function(VectorFunctionSpace(mesh, "DG", 0))
Orientation.assign(as_vector([1.0, 0.0]))

un = dot(Orientation, n)

# define the function space 
V_h = FunctionSpace(mesh, 'DG', degree)
W_h = VectorFunctionSpace(mesh, 'DG', degree)
P_h = TensorFunctionSpace(mesh, 'DG', degree)
R = FunctionSpace(mesh, 'R', 0)

# define the intial value
u_exact_func = as_vector([sin(x)*cos(y)*exp(-2*tc), -cos(x)*sin(y)*exp(-2*tc)]) 
p_exact_func = 0.25*(cos(2*x)+cos(2*y))*exp(-4*tc)

u = Function(W_h)
p = Function(V_h)
u = Function(W_h).interpolate(u_exact_func)
p = Function(V_h).interpolate(p_exact_func)


# Step 1
Mixed_space_1 = W_h * P_h * W_h
u_tilde_trial, Q_trial, r_n_trial = TrialFunctions(Mixed_space_1)
theta_1, Theta, theta_2  = TestFunctions(Mixed_space_1)

Q_flux = conditional(un('+') > 0, Q_trial('-'), Q_trial('+'))
u_tilde_trial_flux = conditional(un('+') > 0, u_tilde_trial('+'), u_tilde_trial('-'))
p_flux = conditional(un('+') > 0, p('-'), p('+'))

# Eq.(1.1)
LHS_1 = dot(u_tilde_trial, theta_1)*dx \
      + dtc*( dot(dot(u, nabla_grad(u_tilde_trial)), theta_1)*dx
            + nu* inner(Q_trial,grad(theta_1))*dx 
            - nu* dot(dot(Q_flux,n('+')), theta_1('+'))*dS 
            - nu* dot(dot(Q_flux,n('-')), theta_1('-'))*dS)\
      + dtc* dot(r_n_trial, theta_1)*dx
RHS_1 = dot(u,theta_1)*dx

# Eq.(1.2)
LHS_1 += inner(Q_trial, Theta)*dx \
       + dot(u_tilde_trial, div(Theta))*dx \
       - dot(u_tilde_trial_flux, dot(Theta('+'), n('+'))) *dS \
       - dot(u_tilde_trial_flux, dot(Theta('-'), n('-'))) *dS 

# Eq.(1.3)
LHS_1 += dot(r_n_trial, theta_2)*dx 
RHS_1 += - p * div(theta_2)*dx \
         + p_flux* jump(theta_2, n)*dS

Solution_Mixed_space_1 = Function(Mixed_space_1)
prob1 = LinearVariationalProblem(LHS_1, RHS_1, Solution_Mixed_space_1)
solv1 = LinearVariationalSolver(prob1)

solv1.solve()

u_tilde = Function(W_h)
r_n = Function(W_h)


u_tilde.assign(Solution_Mixed_space_1.sub(0))
r_n.assign(Solution_Mixed_space_1.sub(2))

# Step 2
Mixed_space_2_1 = W_h * W_h * V_h * R 
Solver_Lagrange_multiplier(Mixed_space_2_1, un, r_n, u_tilde)


Mixed_space_2 = W_h * W_h * V_h 
nullmodes = MixedVectorSpaceBasis(
        Mixed_space_2, [Mixed_space_2.sub(0),VectorSpaceBasis(constant=True),VectorSpaceBasis(constant=True)])
Solver_nullspace(Mixed_space_2, un, r_n, u_tilde, nullmodes)

Output:

R_space, u_max: 1.0000052377552149
R_space, p_max: 0.5000000109023262
nullspace solver, u_max: 1.0147435123510231
nullspace solver, p_max: 1291.2431613978533
No nullspace, u_max: 1.0147435123454347
No nullspace, p_max: 1276.4156866953888

case 2

Dirichlet boundary conditions for both velocity $\boldsymbol{u}$ and pressure $p$, i.e.

$\left.\boldsymbol{u}\right|_{\partial \Omega}=0 $

$\left.p\right|{\partial \Omega}=p{\text{exact}} =\frac{1}{4}e^{-4t}\left(\cos(2x) + \cos(2y)\right)$

Question

This boundary condition should be well-posed, but the pressure $p^{n+1}$ still blows up, and the result is correct after I add the Lagrange multipliers.

case 3

$\boldsymbol{u}$ is periodic boundary condition and pressure $p$ is Dirichlet boundary condition

$\left.p\right|{\partial \Omega}=p{\text{exact}} =\frac{1}{4}e^{-4t}\left(\cos(2x) + \cos(2y)\right)$

Question
Here I do not know how to define the mesh and FunctionSpace.