Basic Simulation Code

OpenCMP includes large amounts of auxilliary code to handle interactions with the user interface and the various post-processing capabilities. As such, it is recommended that new models be developed and debugged as short stand-alone scripts before being incorporated into OpenCMP.

The following code solves the incompressible Navier-Stokes equations using IMEX-style linearization and a first-order IMEX time discretization. It can be used as an example and a starting point for further model development.

# Import the necessary NGSolve/Netgen modules.
import ngsolve as ngs
from netgen.geom2d import SplineGeometry

# Generate the mesh.
geo = SplineGeometry()
geo.AddRectangle((0,0), (2,0.41), bcs=['wall', 'outlet', 'wall', 'inlet'])
geo.AddCircle((0.2,0.2), r=0.05, leftdomain=0, rightdomain=1, bc='cyl')
mesh = ngs.Mesh(geo.GenerateMesh(maxh=0.05))
mesh.Curve(3)

# Define the time variables.
t_0 = 0.0
t_end = 1.0
dt = 0.001

t_param = ngs.Parameter(t_0)
dt_param = ngs.Parameter(dt)

# Various problem values.
kv = 1e-3
interp_ord = 3
u_inlet = ngs.CoefficientFunction((6 * ngs.y * (0.41 - ngs.y) / (0.41 * 0.41), 0.0))
h_outlet = ngs.CoefficientFunction((0.0, 0.0))
f = ngs.CoefficientFunction((0.0, 0.0))
n = ngs.specialcf.normal(mesh.dim)

# Define the finite element space.
V = ngs.VectorH1(mesh, order=interp_ord, dirichlet='wall|inlet|cyl')
Q = ngs.H1(mesh, order=interp_ord - 1)
X = ngs.FESpace([V,Q])

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

# Define gridfunctions to hold the current solution and previous time step.
gfu = ngs.GridFunction(X)
gfu_0 = ngs.GridFunction(X)

# Obtain the initial condition from a Stokes solve.
a_st = ngs.BilinearForm(X)
a_st += (kv * ngs.InnerProduct(ngs.grad(u), ngs.grad(v)) - ngs.div(u) * q - ngs.div(v) * p - 1e-10 * p * q) * ngs.dx

L_st = ngs.LinearForm(X)
L_st += v * f * ngs.dx
L_st += -v * h_outlet * ngs.ds(definedon=mesh.Boundaries('outlet'))

c_st = ngs.Preconditioner(type="direct", bf=a_st, flags = {"inverse" : "umfpack" })
a_st.Assemble()
L_st.Assemble()

gfu.components[0].Set(u_inlet, definedon=mesh.Boundaries('inlet'))
ngs.BVP(bf=a_st, lf=L_st, gf=gfu, pre=c_st, maxsteps=100, prec=1e-10).Do()

gfu_0.vec.data = gfu.vec

# Construct the INS bilinear and linear forms.
a = ngs.BilinearForm(X)
a += (u * v / dt_param) * ngs.dx
a += (kv * ngs.InnerProduct(ngs.grad(u), ngs.grad(v)) - ngs.div(u) * q - ngs.div(v) * p - 1e-10 * p * q) * ngs.dx

L = ngs.LinearForm(X)
L += (gfu_0.components[0] * v / dt_param) * ngs.dx
L += (v * f - ngs.InnerProduct(ngs.grad(gfu_0.components[0]) * gfu_0.components[0], v)) * ngs.dx
L += v * h_outlet * ngs.ds(definedon=mesh.Boundaries('outlet'))

# Construct the INS preconditioner.
c = ngs.Preconditioner(type="direct", bf=a, flags = {"inverse" : "umfpack" })

# Iterate through time.
while t_param.Get() < t_end:
    t_param.Set(t_param.Get() + dt_param.Get())

    # Assemble the bilinear and linear forms and update the preconditioner.
    a.Assemble()
    L.Assemble()
    c.Update()

    # Apply the Dirichlet boundary conditions.
    gfu.components[0].Set(u_inlet, definedon=mesh.Boundaries('inlet'))

    # Single solve.
    ngs.BVP(bf=a, lf=L, gf=gfu, pre=c, maxsteps=20, prec=1e-10).Do()

    # Update the value of the previous time step.
    gfu_0.vec.data = gfu.vec