7. Navier-Stokes on moving domains

Martin Ehrmann, Dominik Freinberger, Felix Strasser

CSE - student seminar SS22, TU Wien

# NGSolve
from ngsolve import *
from netgen.occ import *
from ngsolve.webgui import Draw
from netgen.webgui import Draw as DrawGeo

# Timing
from time import sleep

# Widgets and Sliders
import ipywidgets as widgets
from ipywidgets import interact, fixed

7.1. Geometry and Mesh

# create geometry
square = MoveTo(0,0).Rectangle(5,1).Face()
hole = Circle((2.5,0.5), 0.2).Face()
hole.edges.name = "circle"
hole.edges.maxh = 0.05
square.edges.Max(Y).name = "top"
square.edges.Min(Y).name = "bottom"
square.edges.Min(X).name = "inlet"
square.edges.Max(X).name = "outlet"
shape = square - hole
#DrawGeo(shape)

# generate mesh out of geometry
geo = OCCGeometry(shape, dim=2)
mesh = Mesh(geo.GenerateMesh(maxh=0.1)).Curve(4)
Draw(mesh);

7.2. Wall and Circle Deformation

# define finite element space for deformation function
fes_def = VectorH1(mesh,order=3,dirichlet=".*")

# trial and test functions, bilinear form
u_def, v_def = fes_def.TnT()
a_def = BilinearForm(InnerProduct(grad(u_def), grad(v_def))*dx).Assemble()
A_def_inv = a_def.mat.Inverse(freedofs=fes_def.FreeDofs())

# define grid functinos that later
# will hold the deformations
gf_def_top = GridFunction(fes_def) # upper boundary
gf_def_bot = GridFunction(fes_def) # lower boundary
gf_def_rad = GridFunction(fes_def) # radius of hole
gf_def = GridFunction(fes_def) # this will hold all deformations

# initial values for deformations
delta_top = Parameter(-0.005)
delta_bot = Parameter(0.005)
delta_rad = Parameter(0.02)

# update grid functions with actual deformations 
gf_def_top.Set((0, 2*x*(5-x)), definedon=mesh.Boundaries("top"))
gf_def_bot.Set((0, 2*x*(5-x)), definedon=mesh.Boundaries("bottom"))
gf_def_rad.Set(5*CF((x-2.5, y-0.5)), definedon=mesh.Boundaries("circle"))
gf_def_top.vec.data -= A_def_inv@a_def.mat*gf_def_top.vec
gf_def_bot.vec.data -= A_def_inv@a_def.mat*gf_def_bot.vec
gf_def_rad.vec.data -= A_def_inv@a_def.mat*gf_def_rad.vec
# combine all deformations into one grid function
gf_def.vec.data = delta_top.Get() * gf_def_top.vec + \
                  delta_bot.Get() * gf_def_bot.vec + \
                  delta_rad.Get() * gf_def_rad.vec

7.3. Mean Value Deformation (Bend of the Pipe)

mesh.SetDeformation(None)
# define finite element space for deformation function
fes_def_2 = VectorH1(mesh,order=3,dirichlet="inlet")
fes_number = NumberSpace(mesh)
X_def = fes_def_2 * fes_number

# Parameter here
delta_out = Parameter(0.3)

(u_def_2, lam), (v_def_2, mu) = X_def.TnT()

a_def_2 = BilinearForm(X_def)
a_def_2 += InnerProduct(Sym(grad(u_def_2)), Sym(grad(v_def_2)))*dx
a_def_2 += (u_def_2[1]*mu + v_def_2[1]*lam)*ds("outlet")
a_def_2.Assemble()

A_def_inv_2 = a_def_2.mat.Inverse(freedofs=X_def.FreeDofs())

gf_mean_out = GridFunction(X_def)
l_def_2 = LinearForm(X_def)
l_def_2 += 1*mu*ds("outlet")
l_def_2.Assemble()

gf_mean_out.vec.data = A_def_inv_2*l_def_2.vec 

# add deformation to "global" deformation function
gf_def.vec.data += delta_out.Get()*gf_mean_out.components[0].vec

#Draw(gf_mean_out.components[0])
Draw(gf_def, deformation=gf_def);

7.4. Navier Stokes

mesh.SetDeformation(None)

# timestepping parameters
t = 0; i = 0
tend = 1000
tau = 0.001 # timestep

# define FEM space
order=4
VT = HDiv(mesh, order=order, dirichlet="top|bottom|inlet|circle")
VF = TangentialFacetFESpace(mesh, order=order, dirichlet="top|bottom|inlet|circle|outlet")
Q = L2(mesh, order=order-1)
X = VT*VF*Q

# trial and test functions
u, uhat, p = X.TrialFunction()
v, vhat, q = X.TestFunction()

n = specialcf.normal(mesh.dim)
h = specialcf.mesh_size
dS = dx(element_boundary=True)

# viscosity
nu = 1e-3

def tang(vec):
    return vec - (vec*n)*n

# Thesis Christoph Lehrenfeld, page 71
stokes = nu*InnerProduct(Grad(u), Grad(v)) * dx \
         + nu*InnerProduct(Grad(u)*n, tang(vhat-v)) * dS \
         + nu*InnerProduct(Grad(v)*n, tang(uhat-u)) * dS \
         + nu*4*order*order/h * InnerProduct(tang(vhat-v), tang(uhat-u)) * dS \
         + div(u)*q*dx + div(v)*p*dx -1e-10/nu*p*q*dx

# define and assemble bilinear form
a = BilinearForm(stokes).Assemble()

# assemble linear form
f = LinearForm(X).Assemble()

# grid function for the solution
gfu = GridFunction(X)

# inflow at inlet
uin = CoefficientFunction( (1*1.5*4*y*(1-y), 0) )
gfu.components[0].Set(uin, definedon=mesh.Boundaries("inlet"))

# solve Stokes problem for initial conditions:
inv_stokes = a.mat.Inverse(X.FreeDofs(), inverse="sparsecholesky")

res = f.vec.CreateVector()
res.data = f.vec - a.mat*gfu.vec
gfu.vec.data += inv_stokes * res

# matrix for implicit Euler 
mstar = BilinearForm(X)
mstar += u*v*dx+tau*stokes
mstar.Assemble()
inv = mstar.mat.Inverse(X.FreeDofs(), inverse="sparsecholesky")

# the non-linear term
conv = BilinearForm(X, nonassemble = True)
conv += -InnerProduct(Grad(v)*u, u)*dx # (bonus_intorder=order)

un = u*n
# upwind = IfPos(un, un*u, un*u.Other(uin))
upwind = IfPos(un, un*u, un*u.Other(bnd=mesh.BoundaryCF({"inlet":uin}, default=0.5*u)))
conv += InnerProduct (upwind, v) * dx(element_boundary=True)


VL2 = VectorL2(mesh, order=order, piola=True)
ul2,vl2 = VL2.TnT()
conv_l2 = BilinearForm(VL2, nonassemble=True)
conv_l2 += InnerProduct(grad(vl2)*ul2, ul2).Compile() * dx
conv_l2 += (-IfPos(ul2 * n, ul2*n*ul2*vl2, ul2*n*ul2.Other(bnd=mesh.BoundaryCF({"inlet":uin}, default=0.5*u))*vl2)).Compile() * dS

convertl2 = VT.ConvertL2Operator(VL2) @ X.Restriction(0)
conv_operator = convertl2.T @ conv_l2.mat @ convertl2



# visualization
scene = Draw(gfu.components[0], mesh, order=3) 

# global dictionary holding deformation parameter values
p_dict = {"p1":0, "p2":0, "p3":0, "p4":0}

# function to be executed by interactive() method
# used to update the slider values queried by the Navier-Stokes thread
def parameter_function(p1, p2, p3, p4):
    p_dict["p1"] = p1
    p_dict["p2"] = p2
    p_dict["p3"] = p3
    p_dict["p4"] = p4
    return 0 

# the slider widget
params = widgets.interactive(parameter_function, 
                             p1=widgets.FloatSlider(value=0, min=-0.02, max=0.02, step=0.001, 
                                                    description='upper curve:', disabled=False, 
                                                    continuous_update=True, orientation='horizontal',
                                                    readout=True, readout_format='.3f'),
                             p2=widgets.FloatSlider(value=0, min=-0.02, max=0.02, step=0.001, 
                                                    description='lower curve:', disabled=False, 
                                                    continuous_update=True, orientation='horizontal',
                                                    readout=True, readout_format='.3f'),
                             p3=widgets.FloatSlider(value=0, min=-0.1, max=0.1, step=0.002, 
                                                    description='circle radius:', disabled=False, 
                                                    continuous_update=True, orientation='horizontal',
                                                    readout=True, readout_format='.3f'),
                             p4=widgets.FloatSlider(value=0, min=-1, max=1, step=0.02, 
                                                    description='pipe bend:', disabled=False, 
                                                    continuous_update=True, orientation='horizontal',))

# this is the main function to be run on the separate thread
# it solves the navier stokes problem and repeatedly queries
# the slider values on the main thread (executed only 1x)
def navier_stokes():
    # initial visualization
    scene = Draw(gfu.components[0], mesh, order=3, \
                vectors = { "grid_size":40})
    
    t = 0
    i = 0
    
    # store initial slider values
    p1, p2, p3, p4 = p_dict["p1"], p_dict["p2"], p_dict["p3"], p_dict["p4"]

    # update pipe deformation
    delta_out.Set(p4)
    # update radius deformation
    delta_rad.Set(p3)
    # update upper deformation
    delta_top.Set(p1)
    # update lower deformation
    delta_bot.Set(p2)

    # apply updated parameters to grid function
    gf_def.vec.data = delta_top.Get() * gf_def_top.vec + \
                      delta_bot.Get() * gf_def_bot.vec + \
                      delta_rad.Get() * gf_def_rad.vec + \
                      delta_out.Get() * gf_mean_out.components[0].vec
    
    # apply deformation function to mesh
    mesh.SetDeformation(gf_def)
    
    # assemble RHS and LHS and Euler matrix
    a.Assemble()
    f.Assemble()
    mstar.Assemble()
    inv.Update()
    
    # update convection term
    conv.Apply (gfu.vec, res)
    res.data += a.mat*gfu.vec
    gfu.vec.data -= tau * inv * res   
    
    # update plot
    scene.Redraw()
    
    # the main Navier-Stokes timestepping loop
    SetNumThreads(4)
    with TaskManager(pajetrace=10**8):
        while t < tend:
            # only update parameters and re-assemble system 
            # if slider values have changed
            p1_changed = (p_dict["p1"] != p1)
            p2_changed = (p_dict["p2"] != p2)
            p3_changed = (p_dict["p3"] != p3)
            p4_changed = (p_dict["p4"] != p4)

            if p1_changed or p2_changed or p3_changed or p4_changed:

                # update slider values
                p1, p2, p3, p4 = p_dict["p1"], p_dict["p2"], p_dict["p3"], p_dict["p4"]
                # update pipe deformation
                delta_out.Set(p4)
                # update radius deformation
                delta_rad.Set(p3)
                # update upper deformation
                delta_top.Set(p1)
                # update lower deformation
                delta_bot.Set(p2)

                # apply updated parameters to grid function
                gf_def.vec.data = delta_top.Get() * gf_def_top.vec + \
                                  delta_bot.Get() * gf_def_bot.vec + \
                                  delta_rad.Get() * gf_def_rad.vec + \
                                  delta_out.Get() * gf_mean_out.components[0].vec

                # apply deformation function to mesh
                mesh.SetDeformation(gf_def)

                # re-assemble RHS and LHS and Euler matrix
                a.Assemble()
                f.Assemble()
                mstar.Assemble()
                inv.Update()

            # update convection term
            # conv.Apply (gfu.vec, res)
            # res.data = conv.mat * gfu.vec
            res.data = -conv_operator * gfu.vec
            res.data += a.mat*gfu.vec
            gfu.vec.data -= tau * inv * res   

            # do a time step
            t = t + tau; i = i + 1
            # update plot
            if i%10 == 0: scene.Redraw()
            #if i%500 == 0: print(p_dict)

import threading

# create separate thread with main function "navier_stokes" running the simulation
thread = threading.Thread(target=navier_stokes)

# activate the sliders on the main thread
display(params)

# start the separate thread
thread.start()