TS for Stokes

[JaroslavHron]

Hi, I am trying to use the DAEProblem/DAESolver for the mixed Stokes problem with Taylor-Hood discretization.

Setup is a simple Poiseuille flow with parabolic inflow (one has analytic solution, which is included in the FE space).

However, the DAESolver does not seem to work. Even if I set the analytic solution as the initial data, the initial SNES residuum is large. Am I doing something wrong?

The output of the script:

Solving steady problem of size: W.dim()=2427
  0 SNES Function norm 7.312824288266e+00
  1 SNES Function norm 1.062331577291e-13
Nonlinear firedrake_0_ solve converged due to CONVERGED_FNORM_RELATIVE iterations 1
Steady state solution errors ev=3.300400397290283e-15  ep=1.4578292646079307e-13
Solving timedep. problem of size: W.dim()=2427
Report solution at step=0   time=0.0
Errors are ev=3.300400397290283e-15  ep=1.4578292646079307e-13
0 TS dt 0.01 time 0.
    0 SNES Function norm 2.921282146651e+00
    1 SNES Function norm 1.177247644806e+02
    2 SNES Function norm 2.929534951870e+00
    3 SNES Function norm 3.923617803673e+01
    4 SNES Function norm 2.954809771147e+00
    5 SNES Function norm 2.354003881745e+01
    6 SNES Function norm 2.997444603627e+00
    7 SNES Function norm 1.681310297074e+01
    8 SNES Function norm 3.058025333099e+00
    9 SNES Function norm 1.307622541496e+01

The script:

from firedrake import *
from firedrake.petsc import PETSc
import firedrake_ts

# Data
mu = Constant(1e0)
V_max = Constant(1.0)
L = 2.0
D = 0.5

# Build mesh
mesh = RectangleMesh( 8, 16, D/2.0, L/2.0, originX=-D/2.0, originY=-L/2.0, diagonal='crossed')

# Taylor-Hood finite elements
Ep = FiniteElement("CG",mesh.ufl_cell(),1)
Ev = VectorElement("CG",mesh.ufl_cell(),2)
W = FunctionSpace(mesh,MixedElement([Ev, Ep]))

# Exact solution
x, y = SpatialCoordinate(mesh)
vex = as_vector([0, 4.0*V_max*((D/2.0)*(D/2.0) - x*x)/(D*D)])
pex = -(V_max*8.0*mu/(D*D))*(y - L/2.0)

v_in = interpolate(vex, W.sub(0))

bcs = [
    DirichletBC(W.sub(0), v_in, 3),
    DirichletBC(W.sub(0), Constant((0,0)), 1),
    DirichletBC(W.sub(0), Constant((0,0)), 2)
] 

v_, p_ = TestFunctions(W)
w = Function(W)
v, p = split(w)

F = mu*inner(grad(v), grad(v_))*dx - inner(p, div(v_))*dx + inner(div(v), p_)*dx
J = derivative(F,w)

# First get the steady state solution
print(f"Solving steady problem of size: {W.dim()=}")
problem=NonlinearVariationalProblem(F,w,bcs,J) 
lu = {
    "snes_type": "newtonls",
    "snes_monitor": None,
    "snes_converged_reason": None,
    "ksp_type": "preonly",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": "mumps"
}
solver = NonlinearVariationalSolver(problem, solver_parameters=lu) 
solver.solve()

v, p = w.subfunctions
ev=sqrt(assemble( (v-vex)**2 *dx ) )
ep=sqrt(assemble( (p-pex)**2 *dx ) )
print(f"Steady state solution errors {ev=}  {ep=}")

# Add the time derivative
wdot = Function(W)
vdot, _ = split(w)
F += inner(vdot,v_)*dx

# Try the TS time solver

def monitor(ts, step, time, x):
    PETSc.Sys.Print(f"Report solution at {step=}   {time=}")    
    ctx=ts.getAppCtx()
    w=ctx['w']
    v, p = w.subfunctions
    ev=sqrt(assemble( (v-vex)**2 *dx ) )
    ep=sqrt(assemble( (p-pex)**2 *dx ) )
    print(f"Errors are {ev=}  {ep=}")

lu = {
    "ts_monitor": None,
    "ts_type": 'beuler',
    "ts_max_snes_failures": 1,
    "ts_dt": 0.01,
    "ts_exact_final_time": "stepover",
    "snes_type": "newtonls",
    "snes_monitor": None,
    "snes_converged_reason": None,
    "snes_max_it": 200,
    "snes_rtol": 1e-10,
    "snes_atol": 1e-10,
    "snes_linesearch_type": "cp",
    "ksp_type": "preonly",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": "mumps"
}

problem = firedrake_ts.DAEProblem(F, w, wdot, (0,1), bcs=bcs)
solver = firedrake_ts.DAESolver(problem, solver_parameters=lu, options_prefix='TSNS', monitor_callback=monitor)

solver.ts.setProblemType(solver.ts.ProblemType.NONLINEAR)
solver.ts.setEquationType(solver.ts.EquationType.DAE_IMPLICIT_INDEX2)

ctx={ 'w': w }
solver.ts.setAppCtx(ctx)

print(f"Solving timedep. problem of size: {W.dim()=}")
solver.solve()

[IvanYashchuk]

Hi! It's been a while since I ran a Stokes solver myself. I will take a look at this over the weekend. In the meanwhile, I suggest posting a link to this issue in the Firedrake Slack there are many smart and helpful people.

[JaroslavHron]

Hi, thanks, actually there is a silly typo in the code above: vdot, _ = split(w) should be vdot, _ = split(wdot)... which was part of the problem. Correcting that and with the updates in #20 the stokes problem works.

The last issue which remains is a time dependent DirichletBC - can this be done? Something like this

t = Constant(0.0)
v_in = as_vector([4.0*t, 0.0])
problem = firedrake_ts.DAEProblem(F, w, wdot, (0,10.0), bcs=bcs, time=t)

doesn't seem to work....

[IvanYashchuk]

Nice that you got it working!

One way to get time-dependent BC is to include the condition in the weak form using Nitsche's method. Daniel Shapero has a post about this topic with example Firedrake code: https://shapero.xyz/posts/nitsches-method-nonlinear/

It seems that currently time= constant is unused: https://github.com/IvanYashchuk/firedrake-ts/blob/1186966053ff7ce04b80bc15c02f91f35c0ba476/firedrake_ts/ts_solver.py#L101 I think it should be possible to fetch the current time from PETSc TS (https://petsc.org/release/manualpages/TS/TSGetTime/), update the self.time constant with the fetched value, and do this right before the application of BCs here: https://github.com/IvanYashchuk/firedrake-ts/blob/1186966053ff7ce04b80bc15c02f91f35c0ba476/firedrake_ts/ts_solver.py#L341 Would you try that?