Using Poiseuille example (Test problem 1: Channel flow) with more complex geometries generated with gmsh and additional questions regarding units, implementation of Power Law for viscosity, and best approach for what I'm trying to accomplish

zuehre · January 25, 2024, 9:09pm

I’m trying to set up a simulation using Test problem 1: Channel flow (Poiseuille flow) and gmsh to generate 2d nozzle meshes. Below is the code I am using:

import os
import numpy as np
import matplotlib.pyplot as plt
import gmsh
import pyvista
from lib import dimdict
import tqdm.autonotebook
from mpi4py import MPI
from petsc4py import PETSc
from ufl import (FacetNormal, FiniteElement, Identity, Measure, TestFunction, 
                 TrialFunction, VectorElement, as_vector, div, dot, ds, dx, inner, 
                 lhs, grad, nabla_grad, rhs, sym, sqrt)

from dolfinx.plot import vtk_mesh
from dolfinx.fem import (Constant, Function, FunctionSpace,
                         assemble_scalar, dirichletbc, form, locate_dofs_topological, 
                         set_bc, locate_dofs_geometrical)

from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,
                               create_vector, create_matrix, set_bc)

from dolfinx.graph import adjacencylist
from dolfinx.geometry import bb_tree, compute_collisions_points, compute_colliding_cells
from dolfinx.io import (XDMFFile, VTXWriter, distribute_entity_data, gmshio)
from dolfinx.mesh import create_mesh, meshtags_from_entities
from ufl import (FacetNormal, FiniteElement, MixedElement, Identity, Measure, TestFunction,
                 TestFunctions, TrialFunction, VectorElement, as_vector, div, dot, ds, 
                 dx, inner, lhs, grad, nabla_grad, rhs, sym, split)

## Creating mesh for a syringe
comm = MPI.COMM_WORLD
model_rank = 0
syringe_param_entry = 0
dimdictent=dimdict.ents['Entries'][syringe_param_entry]

# Creating gmsh mesh
gdim = 2
vals = dimdictent["Params"]
gmsh.initialize()
# gmsh.option.setNumber("Mesh.Algorithm", 5)

# Syringe dimensions
syringe_radius = vals["syringe_radius"]["Value"]
syringe_length = vals["syringe_length"]["Value"]

# Nozzle dimensions
nozzle_base_radius = vals["nozzle_base_radius"]["Value"]
nozzle_length = vals["nozzle_length"]["Value"]
nozzle_tip_radius = vals["nozzle_tip_radius"]["Value"]

# Nozzle stage dimensions
nozzle_stage_length = vals["nozzle_stage_length"]["Value"]
nozzle_stage_radius = vals["nozzle_stage_radius"]["Value"]

# Outline points
stagelenx = syringe_length + nozzle_stage_length
nozzlelenx = stagelenx + nozzle_length

# Define mesh sizes
ms1 = syringe_radius/2
ms2 = nozzle_base_radius/2
ms3 = nozzle_tip_radius
gmsh.option.setNumber("Mesh.MeshSizeMax", 30)

# Order points
if nozzle_stage_length > 0:
    ptlst = [
        (0,syringe_radius),
        (syringe_length,syringe_radius),
        (syringe_length,nozzle_stage_radius),
        (stagelenx,nozzle_stage_radius),
        (stagelenx,nozzle_base_radius),
        (nozzlelenx,nozzle_tip_radius),
        (nozzlelenx,-nozzle_tip_radius),
        (stagelenx,-nozzle_base_radius),
        (stagelenx,-nozzle_stage_radius),
        (syringe_length,-nozzle_stage_radius),
        (syringe_length,-syringe_radius),
        (0,-syringe_radius)
    ]
    nozzle_stage = True
else:
    ptlst = [
        (0,syringe_radius),
        (syringe_length,syringe_radius),
        (syringe_length,nozzle_base_radius),
        (nozzlelenx,nozzle_tip_radius),
        (nozzlelenx,-nozzle_tip_radius),
        (syringe_length,-nozzle_base_radius),
        (syringe_length,-syringe_radius),
        (0,-syringe_radius)
    ]
    nozzle_stage = False

# Defining starting point
gmsh.model.occ.addPoint(0,syringe_radius,0,ms1,1)

# Defining gmsh points from precompiled list and assigning resolution based on radius
count = 1
lnlst = []
for pt in range(len(ptlst)-1):
    
    if abs(ptlst[pt+1][1]) == syringe_radius:
        res = ms1
    elif abs(ptlst[pt+1][1]) == nozzle_stage_radius:
        res = ms2
    elif abs(ptlst[pt+1][1]) == nozzle_base_radius:
        res = ms2
    elif abs(ptlst[pt+1][1]) == nozzle_tip_radius:
        res = ms3

    if ptlst[pt] != ptlst[pt+1]:
        gmsh.model.occ.addPoint(ptlst[pt+1][0],ptlst[pt+1][1],0,res,count+1)
        gmsh.model.occ.addLine(count,count+1,count)
        lnlst.append(count)
        count += 1

# Defining ending point and last connecting line
gmsh.model.occ.addPoint(0,-syringe_radius,0,ms1,count+1)
gmsh.model.occ.addLine(count,1,count)
lnlst.append(count)

# Creating surface
gmsh.model.occ.addCurveLoop(lnlst,lnlst[-1]+1) # base surface
gmsh.model.occ.addPlaneSurface([lnlst[-1]+1],lnlst[-1]+2)
    
# Synchronizing mesh parts
gmsh.model.occ.synchronize()

# Defining nodes where fluid will flow
surfaces = gmsh.model.getEntities(dim=gdim)
assert (len(surfaces) == 1)
fm = surfaces[0][1]
gmsh.model.addPhysicalGroup(surfaces[0][0], [fm], fm)
gmsh.model.setPhysicalName(surfaces[0][0], fm, "Fluid")

# Finding boundary, inflow, and outflow nodes
inlet_marker, outlet_marker, wall_marker, idk_marker = fm+1,fm+2,fm+3,fm+4
inflow, outflow, walls, leftovers = [], [], [], []
boundaries = gmsh.model.getBoundary(surfaces,oriented=False)
for boundary in boundaries:
    center_of_mass = gmsh.model.occ.getCenterOfMass(boundary[0], boundary[1])
    if np.allclose(center_of_mass, [0, 0, 0]):
        inflow.append(boundary[1])
    elif np.allclose(center_of_mass, [nozzlelenx, 0, 0]):
        outflow.append(boundary[1])
    else:
        walls.append(boundary[1])

gmsh.model.addPhysicalGroup(1, walls, wall_marker)
gmsh.model.setPhysicalName(1, wall_marker, "Walls")
gmsh.model.addPhysicalGroup(1, inflow, inlet_marker)
gmsh.model.setPhysicalName(1, inlet_marker, "Inlet")
gmsh.model.addPhysicalGroup(1, outflow, outlet_marker)
gmsh.model.setPhysicalName(1, outlet_marker, "Outlet")
gmsh.model.addPhysicalGroup(1, leftovers, idk_marker)
gmsh.model.setPhysicalName(1, idk_marker, "IDK")

# Generating the mesh
gmsh_instance = gmsh.model.mesh.generate(gdim)

# Creating mesh model with gmsh mesh
mesh, ct, ft = gmshio.model_to_mesh(gmsh.model,comm=comm,rank=comm.rank,gdim=gdim)
ft.name = "Facet Markers"
ct.name = "Cell Markers"

with XDMFFile(mesh.comm, "results/mesh.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)

# Finalize gmsh
gmsh.finalize()

# The following resources were used to build this script:
# https://jsdokken.com/dolfinx-tutorial/chapter2/ns_code2.html
# https://fenicsproject.org/qa/11271/nonlinear-non-newtonian-navier-stokes/
# https://fenicsproject.discourse.group/t/ipcs-solver-with-time-dependent-pressure-bc/10616/5
# https://jsdokken.com/dolfinx-tutorial/chapter2/ns_code1.html

## Defining problem specific parameters
t = 0                       # Start time
T = 10                      # Final time
dt = 1 / 100                # Time step size
num_steps = int(T / dt)

# Physical properties of water (for testing)
# Units are kg, mm, s, mPa, (mm/s), 
RHO = 1     # kg/(mm**3)
MU = 1         # kg/(mm*s)

f = Constant(mesh, PETSc.ScalarType((0, 0)))
k = Constant(mesh, PETSc.ScalarType(dt))
rho = Constant(mesh, PETSc.ScalarType(RHO))     # Density
mu = Constant(mesh, PETSc.ScalarType(MU))      # Viscosity

# Keeping viscosity constant for initial testing, but would like to make it vary with shear rate according to Power Law
# Defining visocity for a shear-thinning fluid using the Power Law
K = Constant(mesh, PETSc.ScalarType(1.0))  # Rheological constant
m = Constant(mesh, PETSc.ScalarType(0.5))  # Power-law index

# Defining inlet / delta pressure (outlet pressure is 0 gauge)
dp = 8           # mPa

# Defining Taylor-Hood function space
v_cg2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
s_cg1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, v_cg2)
Q = FunctionSpace(mesh, s_cg1)

u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)
u_n = Function(V)
u_n.name = "u_n"
U = 0.5 * (u_n + u)
n = FacetNormal(mesh)
fdim = mesh.topology.dim-1

# Define boundary conditions
# Walls
wall_dofs = locate_dofs_topological(V, fdim, ft.find(wall_marker))
u_nonslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)
bcu_walls = dirichletbc(u_nonslip, wall_dofs, V)
bcu = [bcu_walls]

# Inlet
p_inlet = PETSc.ScalarType(dp)
inflow_dofs = locate_dofs_topological(Q, fdim, ft.find(inlet_marker))
bcp_inflow = dirichletbc(p_inlet, inflow_dofs, Q)

# Outlet
p_outlet = PETSc.ScalarType(0)
outflow_dofs = locate_dofs_topological(Q, fdim, ft.find(outlet_marker))
bcp_outlet = dirichletbc(p_outlet, outflow_dofs, Q)
bcp = [bcp_inflow, bcp_outlet]

# Define strain-rate tensor
def epsilon(u):
    return sym(nabla_grad(u))

# Define stress tensor
def sigma(u, p):
    return 2 * mu * epsilon(u) - p * Identity(len(u))

# Define the variational problem for the first step
p_n = Function(Q)
p_n.name = "p_n"
F1 = rho * dot((u - u_n) / k, v) * dx
F1 += rho * dot(dot(u_n, nabla_grad(u_n)), v) * dx
F1 += inner(sigma(U, p_n), epsilon(v)) * dx
F1 += dot(p_n * n, v) * ds - dot(mu * nabla_grad(U) * n, v) * ds
F1 -= dot(f, v) * dx
a1 = form(lhs(F1))
L1 = form(rhs(F1))
A1 = assemble_matrix(a1, bcs=bcu)
A1.assemble()
b1 = create_vector(L1)

# Define variational problem for the second step
u_ = Function(V)
a2 = form(dot(nabla_grad(p), nabla_grad(q)) * dx)
L2 = form(dot(nabla_grad(p_n), nabla_grad(q)) * dx - (rho / k) * div(u_) * q * dx)
A2 = assemble_matrix(a2, bcs=bcp)
A2.assemble()
b2 = create_vector(L2)

# Define variational problem for the third step
p_ = Function(Q)
a3 = form(rho * dot(u, v) * dx)
L3 = form(rho * dot(u_, v) * dx - k * dot(nabla_grad(p_ - p_n), v) * dx)
A3 = assemble_matrix(a3)
A3.assemble()
b3 = create_vector(L3)

# Solver for step 1
solver1 = PETSc.KSP().create(mesh.comm)
solver1.setOperators(A1)
solver1.setType(PETSc.KSP.Type.BCGS)
pc1 = solver1.getPC()
pc1.setType(PETSc.PC.Type.HYPRE)
pc1.setHYPREType("boomeramg")

# Solver for step 2
solver2 = PETSc.KSP().create(mesh.comm)
solver2.setOperators(A2)
solver2.setType(PETSc.KSP.Type.BCGS)
pc2 = solver2.getPC()
pc2.setType(PETSc.PC.Type.HYPRE)
pc2.setHYPREType("boomeramg")

# Solver for step 3
solver3 = PETSc.KSP().create(mesh.comm)
solver3.setOperators(A3)
solver3.setType(PETSc.KSP.Type.CG)
pc3 = solver3.getPC()
pc3.setType(PETSc.PC.Type.SOR)

from pathlib import Path
folder = Path("results")
folder.mkdir(exist_ok=True, parents=True)
vtx_u = VTXWriter(mesh.comm, "dfg2D-3-u.bp", [u_], engine="BP4")
vtx_p = VTXWriter(mesh.comm, "dfg2D-3-p.bp", [p_], engine="BP4")
vtx_u.write(t)
vtx_p.write(t)

## Solving
def u_exact(x):
    values = np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType)
    values[0] = 4 * x[1] * (1.0 - x[1])
    return values

u_ex = Function(V)
u_ex.interpolate(u_exact)
L2_error = form(dot(u_ - u_ex, u_ - u_ex) * dx)

for i in range(num_steps):
    # Update current time step
    t += dt

    # Step 1: Tentative veolcity step
    with b1.localForm() as loc_1:
        loc_1.set(0)
    assemble_vector(b1, L1)
    apply_lifting(b1, [a1], [bcu])
    b1.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b1, bcu)
    solver1.solve(b1, u_.vector)
    u_.x.scatter_forward()

    # Step 2: Pressure corrrection step
    with b2.localForm() as loc_2:
        loc_2.set(0)
    assemble_vector(b2, L2)
    apply_lifting(b2, [a2], [bcp])
    b2.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b2, bcp)
    solver2.solve(b2, p_.vector)
    p_.x.scatter_forward()

    # Step 3: Velocity correction step
    with b3.localForm() as loc_3:
        loc_3.set(0)
    assemble_vector(b3, L3)
    b3.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    solver3.solve(b3, u_.vector)
    u_.x.scatter_forward()
    # Update variable with solution form this time step
    u_n.x.array[:] = u_.x.array[:]
    p_n.x.array[:] = p_.x.array[:]

    # Write solutions to file
    vtx_u.write(t)
    vtx_p.write(t)

    # Compute error at current time-step
    error_L2 = np.sqrt(mesh.comm.allreduce(assemble_scalar(L2_error), op=MPI.SUM))
    error_max = mesh.comm.allreduce(np.max(u_.vector.array - u_ex.vector.array), op=MPI.MAX)
    # Print error only every 20th step and at the last step
    if (i % 20 == 0) or (i == num_steps - 1):
        print(f"Time {t:.2f}, L2-error {error_L2:.2e}, Max error {error_max:.2e}")

# Close xmdf file
vtx_u.close()
vtx_p.close()
b1.destroy()
b2.destroy()
b3.destroy()
solver1.destroy()
solver2.destroy()
solver3.destroy()

def plotfield():
    pyvista.start_xvfb()
    topology, cell_types, geometry = vtk_mesh(V)
    values = np.zeros((geometry.shape[0], 3), dtype=np.float64)
    values[:, :len(u_n)] = u_n.x.array.real.reshape((geometry.shape[0], len(u_n)))

    # Create a point cloud of glyphs
    function_grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
    function_grid["u"] = values
    glyphs = function_grid.glyph(orient="u", factor=0.2)

    # Create a pyvista-grid for the mesh
    grid = pyvista.UnstructuredGrid(*vtk_mesh(mesh, mesh.topology.dim))

    # Create plotter
    plotter = pyvista.Plotter(off_screen=True)
    plotter.add_mesh(grid, style="wireframe", color="k")
    plotter.add_mesh(glyphs)
    plotter.view_xy()
    fig_as_array = plotter.screenshot("glyphs.png")
    return

plotfield()

The dictionary dimdict with geometry information looks like this:

ents=\
{"Entries": 
    [
        {   "Entry ID": 0, "Name": "Test1", "Params": 
            {
                "syringe_radius": {"Value": 60, "Units": "mm"},
                "syringe_length": {"Value": 550, "Units": "mm"},
                "nozzle_base_radius": {"Value": 20, "Units": "mm"},
                "nozzle_length": {"Value": 320, "Units": "mm"},
                "nozzle_tip_radius": {"Value": 0.1, "Units": "mm"},
                "nozzle_stage_length": {"Value": 50, "Units": "mm"},
                "nozzle_stage_radius": {"Value": 30, "Units": "mm"},
            }
        },
        {    "Entry ID": 1, "Name": "Test2", "Params": 
            {
                "syringe_radius": {"Value": 30, "Units": "mm"},
                "syringe_length": {"Value": 50, "Units": "mm"},
                "nozzle_base_radius": {"Value": 20, "Units": "mm"},
                "nozzle_length": {"Value": 20, "Units": "mm"},
                "nozzle_tip_radius": {"Value": 0.1, "Units": "mm"},
                "nozzle_stage_length": {"Value": 10, "Units": "mm"},
                "nozzle_stage_radius": {"Value": 20, "Units": "mm"},
            }
        },
        {    "Entry ID": 2, "Name": "Test3", "Params": 
            {
                "syringe_radius": {"Value": 30, "Units": "mm"},
                "syringe_length": {"Value": 200, "Units": "mm"},
                "nozzle_base_radius": {"Value": 30, "Units": "mm"},
                "nozzle_length": {"Value": 50, "Units": "mm"},
                "nozzle_tip_radius": {"Value": 30, "Units": "mm"},
                "nozzle_stage_length": {"Value": 0, "Units": "mm"},
                "nozzle_stage_radius": {"Value": 0, "Units": "mm"},
            }
        },
        {    "Entry ID": 3, "Name": "Test4", "Params": 
            {
                "syringe_radius": {"Value": 0.5, "Units": "mm"},
                "syringe_length": {"Value": 0.5, "Units": "mm"},
                "nozzle_base_radius": {"Value": 0.5, "Units": "mm"},
                "nozzle_length": {"Value": 0.5, "Units": "mm"},
                "nozzle_tip_radius": {"Value": 0.5, "Units": "mm"},
                "nozzle_stage_length": {"Value": 0, "Units": "mm"},
                "nozzle_stage_radius": {"Value": 0, "Units": "mm"},
            }
        }
    ]
}

When I recreate a unit square mesh (Entry ID 3) and use the same parameters as in the example (mu=1,rho=1,deltapressure=8), the generated plot is comparable to output from the example, though the calculated error is much greater:

Info    : Meshing 1D...
Info    : [  0%] Meshing curve 1 (Line)
Info    : [ 20%] Meshing curve 2 (Line)
Info    : [ 40%] Meshing curve 3 (Line)
Info    : [ 50%] Meshing curve 4 (Line)
Info    : [ 70%] Meshing curve 5 (Line)
Info    : [ 90%] Meshing curve 6 (Line)
Info    : Done meshing 1D (Wall 0.000524887s, CPU 0.000246s)
Info    : Meshing 2D...
Info    : Meshing surface 8 (Plane, Frontal-Delaunay)
Info    : Done meshing 2D (Wall 0.00107404s, CPU 0.001073s)
Info    : 32 nodes 67 elements
Time 0.01, L2-error 1.26e+00, Max error 3.08e+00
Time 0.21, L2-error 1.48e+00, Max error 3.00e+00
Time 0.41, L2-error 1.52e+00, Max error 3.00e+00
Time 0.61, L2-error 1.53e+00, Max error 3.00e+00
Time 0.81, L2-error 1.53e+00, Max error 3.00e+00
Time 1.01, L2-error 1.53e+00, Max error 3.00e+00
Time 1.21, L2-error 1.53e+00, Max error 3.00e+00
Time 1.41, L2-error 1.53e+00, Max error 3.00e+00
Time 1.61, L2-error 1.53e+00, Max error 3.00e+00
Time 1.81, L2-error 1.53e+00, Max error 3.00e+00
Time 2.01, L2-error 1.53e+00, Max error 3.00e+00
Time 2.21, L2-error 1.53e+00, Max error 3.00e+00
Time 2.41, L2-error 1.53e+00, Max error 3.00e+00
Time 2.61, L2-error 1.53e+00, Max error 3.00e+00
Time 2.81, L2-error 1.53e+00, Max error 3.00e+00
Time 3.01, L2-error 1.53e+00, Max error 3.00e+00
Time 3.21, L2-error 1.53e+00, Max error 3.00e+00
Time 3.41, L2-error 1.53e+00, Max error 3.00e+00
Time 3.61, L2-error 1.53e+00, Max error 3.00e+00
Time 3.81, L2-error 1.53e+00, Max error 3.00e+00
Time 4.01, L2-error 1.53e+00, Max error 3.00e+00
Time 4.21, L2-error 1.53e+00, Max error 3.00e+00
Time 4.41, L2-error 1.53e+00, Max error 3.00e+00
Time 4.61, L2-error 1.53e+00, Max error 3.00e+00
Time 4.81, L2-error 1.53e+00, Max error 3.00e+00
Time 5.01, L2-error 1.53e+00, Max error 3.00e+00
Time 5.21, L2-error 1.53e+00, Max error 3.00e+00
Time 5.41, L2-error 1.53e+00, Max error 3.00e+00
Time 5.61, L2-error 1.53e+00, Max error 3.00e+00
Time 5.81, L2-error 1.53e+00, Max error 3.00e+00
Time 6.01, L2-error 1.53e+00, Max error 3.00e+00
Time 6.21, L2-error 1.53e+00, Max error 3.00e+00
Time 6.41, L2-error 1.53e+00, Max error 3.00e+00
Time 6.61, L2-error 1.53e+00, Max error 3.00e+00
Time 6.81, L2-error 1.53e+00, Max error 3.00e+00
Time 7.01, L2-error 1.53e+00, Max error 3.00e+00
Time 7.21, L2-error 1.53e+00, Max error 3.00e+00
Time 7.41, L2-error 1.53e+00, Max error 3.00e+00
Time 7.61, L2-error 1.53e+00, Max error 3.00e+00
Time 7.81, L2-error 1.53e+00, Max error 3.00e+00
Time 8.01, L2-error 1.53e+00, Max error 3.00e+00
Time 8.21, L2-error 1.53e+00, Max error 3.00e+00
Time 8.41, L2-error 1.53e+00, Max error 3.00e+00
Time 8.61, L2-error 1.53e+00, Max error 3.00e+00
Time 8.81, L2-error 1.53e+00, Max error 3.00e+00
Time 9.01, L2-error 1.53e+00, Max error 3.00e+00
Time 9.21, L2-error 1.53e+00, Max error 3.00e+00
Time 9.41, L2-error 1.53e+00, Max error 3.00e+00
Time 9.61, L2-error 1.53e+00, Max error 3.00e+00
Time 9.81, L2-error 1.53e+00, Max error 3.00e+00
Time 10.00, L2-error 1.53e+00, Max error 3.00e+00

So, my first question is what exactly does the error mean and why does it change from 10**-6 to 10**0 with a slightly different mesh geometry?
My second question is, what units need to be used here? In order to achieve dimensionlessness in the Reynolds Number and considering that the mesh is measured out in mm, expressing density as kg/(mm3), viscosity as kg/(mms), and pressure as kg/(mms2) (equivalent to mPa) should work. When using water at 20C as an example, rho would be 1/1000, mu would be 1/1000000, and a reasonable pressure would be 1000 mPa, but these values shoot the error up to 10**93 for the unit square mesh by the second time step:

Info    : Meshing 1D...
Info    : [  0%] Meshing curve 1 (Line)
Info    : [ 20%] Meshing curve 2 (Line)
Info    : [ 40%] Meshing curve 3 (Line)
Info    : [ 50%] Meshing curve 4 (Line)
Info    : [ 70%] Meshing curve 5 (Line)
Info    : [ 90%] Meshing curve 6 (Line)
Info    : Done meshing 1D (Wall 0.000373357s, CPU 0.000192s)
Info    : Meshing 2D...
Info    : Meshing surface 8 (Plane, Frontal-Delaunay)
Info    : Done meshing 2D (Wall 0.0012464s, CPU 0.001246s)
Info    : 32 nodes 67 elements
Time 0.01, L2-error 1.00e+04, Max error 1.00e+04
Time 0.21, L2-error 1.98e+93, Max error 4.94e+93

Similarly, I am also running into stability issues with the simulation when I change the parameters. Specifically, if I increase the linear dimensions of the mesh or use a more complex mesh (not a square), error increases quite a bit. Using a square with linear dimensions increased 10-fold results in the following:

Info    : Meshing 1D...
Info    : [  0%] Meshing curve 1 (Line)
Info    : [ 20%] Meshing curve 2 (Line)
Info    : [ 40%] Meshing curve 3 (Line)
Info    : [ 50%] Meshing curve 4 (Line)
Info    : [ 70%] Meshing curve 5 (Line)
Info    : [ 90%] Meshing curve 6 (Line)
Info    : Done meshing 1D (Wall 0.000636288s, CPU 0.00064s)
Info    : Meshing 2D...
Info    : Meshing surface 8 (Plane, Frontal-Delaunay)
Info    : Done meshing 2D (Wall 0.00112193s, CPU 0.00112s)
Info    : 32 nodes 67 elements
Time 0.01, L2-error 4.47e+05, Max error 1.02e+04
Time 0.21, L2-error 4.47e+05, Max error 1.02e+04
Time 0.41, L2-error 4.47e+05, Max error 1.02e+04
Time 0.61, L2-error 4.47e+05, Max error 1.02e+04
Time 0.81, L2-error 4.47e+05, Max error 1.02e+04
Time 1.01, L2-error 4.47e+05, Max error 1.02e+04
Time 1.21, L2-error 4.47e+05, Max error 1.02e+04
Time 1.41, L2-error 4.47e+05, Max error 1.02e+04
Time 1.61, L2-error 4.47e+05, Max error 1.02e+04
Time 1.81, L2-error 4.47e+05, Max error 1.02e+04
Time 2.01, L2-error 4.47e+05, Max error 1.02e+04
Time 2.21, L2-error 4.47e+05, Max error 1.02e+04
Time 2.41, L2-error 4.47e+05, Max error 1.02e+04
Time 2.61, L2-error 4.47e+05, Max error 1.02e+04
Time 2.81, L2-error 4.47e+05, Max error 1.02e+04
Time 3.01, L2-error 4.47e+05, Max error 1.02e+04
Time 3.21, L2-error 4.47e+05, Max error 1.02e+04
Time 3.41, L2-error 4.47e+05, Max error 1.02e+04
Time 3.61, L2-error 4.47e+05, Max error 1.02e+04
Time 3.81, L2-error 4.47e+05, Max error 1.02e+04
Time 4.01, L2-error 4.47e+05, Max error 1.02e+04
Time 4.21, L2-error 4.47e+05, Max error 1.02e+04
Time 4.41, L2-error 4.47e+05, Max error 1.02e+04
Time 4.61, L2-error 4.47e+05, Max error 1.02e+04
Time 4.81, L2-error 4.47e+05, Max error 1.02e+04
Time 5.01, L2-error 4.47e+05, Max error 1.02e+04
Time 5.21, L2-error 4.47e+05, Max error 1.02e+04
Time 5.41, L2-error 4.47e+05, Max error 1.02e+04
Time 5.61, L2-error 4.47e+05, Max error 1.02e+04
Time 5.81, L2-error 4.47e+05, Max error 1.02e+04
Time 6.01, L2-error 4.47e+05, Max error 1.02e+04
Time 6.21, L2-error 4.47e+05, Max error 1.02e+04
Time 6.41, L2-error 4.47e+05, Max error 1.02e+04
Time 6.61, L2-error 4.47e+05, Max error 1.02e+04
Time 6.81, L2-error 4.47e+05, Max error 1.02e+04
Time 7.01, L2-error 4.47e+05, Max error 1.02e+04
Time 7.21, L2-error 4.47e+05, Max error 1.02e+04
Time 7.41, L2-error 4.47e+05, Max error 1.02e+04
Time 7.61, L2-error 4.47e+05, Max error 1.02e+04
Time 7.81, L2-error 4.47e+05, Max error 1.02e+04
Time 8.01, L2-error 4.47e+05, Max error 1.02e+04
Time 8.21, L2-error 4.47e+05, Max error 1.02e+04
Time 8.41, L2-error 4.47e+05, Max error 1.02e+04
Time 8.61, L2-error 4.47e+05, Max error 1.02e+04
Time 8.81, L2-error 4.47e+05, Max error 1.02e+04
Time 9.01, L2-error 4.47e+05, Max error 1.02e+04
Time 9.21, L2-error 4.47e+05, Max error 1.02e+04
Time 9.41, L2-error 4.47e+05, Max error 1.02e+04
Time 9.61, L2-error 4.47e+05, Max error 1.02e+04
Time 9.81, L2-error 4.47e+05, Max error 1.02e+04
Time 10.00, L2-error 4.47e+05, Max error 1.02e+04

I apologize for any dumb questions; I am not well versed in finite element analysis or how to set the solvers up. Any and all input, specifically suggestions on how to better set up this simulation (modeling shear stress / velocity of pressure-driven flow in complex 2d geometries built with gmsh and using dolfinx to set up the problem and solve) is greatly appreciated.

I am also interested making the viscosity parameter mu update based on the velocity (shear thinning behavior using the Power Law relationship between shear stress and shear rate).

dokken · January 26, 2024, 2:15am

There are many good questions here, that i do not have bandwidth to answer atm.
First of all I would go with a more advanced splitting scheme, as shown in
https://jsdokken.com/dolfinx-tutorial/chapter2/ns_code2.html

Secondly, your mesh is alot coarser. The error of any Finite element simulation is a function of the spatial and temporal representation (dx, dt). How do the error change if you reduce the mesh size?

zuehre · January 26, 2024, 5:12pm

The same error is returned for the unit square when running ~3x more nodes:

Info    : Meshing 1D...
Info    : [  0%] Meshing curve 1 (Line)
Info    : [ 20%] Meshing curve 2 (Line)
Info    : [ 40%] Meshing curve 3 (Line)
Info    : [ 50%] Meshing curve 4 (Line)
Info    : [ 70%] Meshing curve 5 (Line)
Info    : [ 90%] Meshing curve 6 (Line)
Info    : Done meshing 1D (Wall 0.000246468s, CPU 0.000112s)
Info    : Meshing 2D...
Info    : Meshing surface 8 (Plane, Frontal-Delaunay)
Info    : Done meshing 2D (Wall 0.00286078s, CPU 0.002856s)
Info    : 99 nodes 201 elements
Time 0.01, L2-error 1.26e+00, Max error 3.08e+00
Time 0.21, L2-error 1.48e+00, Max error 3.00e+00
Time 0.41, L2-error 1.52e+00, Max error 3.00e+00
Time 0.61, L2-error 1.53e+00, Max error 3.00e+00
Time 0.81, L2-error 1.53e+00, Max error 3.00e+00
Time 1.01, L2-error 1.53e+00, Max error 3.00e+00
Time 1.21, L2-error 1.53e+00, Max error 3.00e+00
Time 1.41, L2-error 1.53e+00, Max error 3.00e+00
Time 1.61, L2-error 1.53e+00, Max error 3.00e+00
Time 1.81, L2-error 1.53e+00, Max error 3.00e+00
Time 2.01, L2-error 1.53e+00, Max error 3.00e+00
Time 2.21, L2-error 1.53e+00, Max error 3.00e+00
Time 2.41, L2-error 1.53e+00, Max error 3.00e+00
Time 2.61, L2-error 1.53e+00, Max error 3.00e+00
Time 2.81, L2-error 1.53e+00, Max error 3.00e+00
Time 3.01, L2-error 1.53e+00, Max error 3.00e+00
Time 3.21, L2-error 1.53e+00, Max error 3.00e+00
Time 3.41, L2-error 1.53e+00, Max error 3.00e+00
Time 3.61, L2-error 1.53e+00, Max error 3.00e+00
Time 3.81, L2-error 1.53e+00, Max error 3.00e+00
Time 4.01, L2-error 1.53e+00, Max error 3.00e+00
Time 4.21, L2-error 1.53e+00, Max error 3.00e+00
Time 4.41, L2-error 1.53e+00, Max error 3.00e+00
Time 4.61, L2-error 1.53e+00, Max error 3.00e+00
Time 4.81, L2-error 1.53e+00, Max error 3.00e+00
Time 5.01, L2-error 1.53e+00, Max error 3.00e+00
Time 5.21, L2-error 1.53e+00, Max error 3.00e+00
Time 5.41, L2-error 1.53e+00, Max error 3.00e+00
Time 5.61, L2-error 1.53e+00, Max error 3.00e+00
Time 5.81, L2-error 1.53e+00, Max error 3.00e+00
Time 6.01, L2-error 1.53e+00, Max error 3.00e+00
Time 6.21, L2-error 1.53e+00, Max error 3.00e+00
Time 6.41, L2-error 1.53e+00, Max error 3.00e+00
Time 6.61, L2-error 1.53e+00, Max error 3.00e+00
Time 6.81, L2-error 1.53e+00, Max error 3.00e+00
Time 7.01, L2-error 1.53e+00, Max error 3.00e+00
Time 7.21, L2-error 1.53e+00, Max error 3.00e+00
Time 7.41, L2-error 1.53e+00, Max error 3.00e+00
Time 7.61, L2-error 1.53e+00, Max error 3.00e+00
Time 7.81, L2-error 1.53e+00, Max error 3.00e+00
Time 8.01, L2-error 1.53e+00, Max error 3.00e+00
Time 8.21, L2-error 1.53e+00, Max error 3.00e+00
Time 8.41, L2-error 1.53e+00, Max error 3.00e+00
Time 8.61, L2-error 1.53e+00, Max error 3.00e+00
Time 8.81, L2-error 1.53e+00, Max error 3.00e+00
Time 9.01, L2-error 1.53e+00, Max error 3.00e+00
Time 9.21, L2-error 1.53e+00, Max error 3.00e+00
Time 9.41, L2-error 1.53e+00, Max error 3.00e+00
Time 9.61, L2-error 1.53e+00, Max error 3.00e+00
Time 9.81, L2-error 1.53e+00, Max error 3.00e+00
Time 10.00, L2-error 1.53e+00, Max error 3.00e+00

And the same case for the unit square when running ~16x more nodes:

Info    : Meshing 1D...
Info    : [  0%] Meshing curve 1 (Line)
Info    : [ 20%] Meshing curve 2 (Line)
Info    : [ 40%] Meshing curve 3 (Line)
Info    : [ 50%] Meshing curve 4 (Line)
Info    : [ 70%] Meshing curve 5 (Line)
Info    : [ 90%] Meshing curve 6 (Line)
Info    : Done meshing 1D (Wall 0.000723788s, CPU 0.000343s)
Info    : Meshing 2D...
Info    : Meshing surface 8 (Plane, Frontal-Delaunay)
Info    : Done meshing 2D (Wall 0.0129035s, CPU 0.012248s)
Info    : 512 nodes 1027 elements
Time 0.01, L2-error 1.26e+00, Max error 3.08e+00
Time 0.21, L2-error 1.48e+00, Max error 3.00e+00
Time 0.41, L2-error 1.52e+00, Max error 3.00e+00
Time 0.61, L2-error 1.53e+00, Max error 3.00e+00
Time 0.81, L2-error 1.53e+00, Max error 3.00e+00
Time 1.01, L2-error 1.53e+00, Max error 3.00e+00
Time 1.21, L2-error 1.53e+00, Max error 3.00e+00
Time 1.41, L2-error 1.53e+00, Max error 3.00e+00
Time 1.61, L2-error 1.53e+00, Max error 3.00e+00
Time 1.81, L2-error 1.53e+00, Max error 3.00e+00
Time 2.01, L2-error 1.53e+00, Max error 3.00e+00
Time 2.21, L2-error 1.53e+00, Max error 3.00e+00
Time 2.41, L2-error 1.53e+00, Max error 3.00e+00
Time 2.61, L2-error 1.53e+00, Max error 3.00e+00
Time 2.81, L2-error 1.53e+00, Max error 3.00e+00
Time 3.01, L2-error 1.53e+00, Max error 3.00e+00
Time 3.21, L2-error 1.53e+00, Max error 3.00e+00
Time 3.41, L2-error 1.53e+00, Max error 3.00e+00
Time 3.61, L2-error 1.53e+00, Max error 3.00e+00
Time 3.81, L2-error 1.53e+00, Max error 3.00e+00
Time 4.01, L2-error 1.53e+00, Max error 3.00e+00
Time 4.21, L2-error 1.53e+00, Max error 3.00e+00
Time 4.41, L2-error 1.53e+00, Max error 3.00e+00
Time 4.61, L2-error 1.53e+00, Max error 3.00e+00
Time 4.81, L2-error 1.53e+00, Max error 3.00e+00
Time 5.01, L2-error 1.53e+00, Max error 3.00e+00
Time 5.21, L2-error 1.53e+00, Max error 3.00e+00
Time 5.41, L2-error 1.53e+00, Max error 3.00e+00
Time 5.61, L2-error 1.53e+00, Max error 3.00e+00
Time 5.81, L2-error 1.53e+00, Max error 3.00e+00
Time 6.01, L2-error 1.53e+00, Max error 3.00e+00
Time 6.21, L2-error 1.53e+00, Max error 3.00e+00
Time 6.41, L2-error 1.53e+00, Max error 3.00e+00
Time 6.61, L2-error 1.53e+00, Max error 3.00e+00
Time 6.81, L2-error 1.53e+00, Max error 3.00e+00
Time 7.01, L2-error 1.53e+00, Max error 3.00e+00
Time 7.21, L2-error 1.53e+00, Max error 3.00e+00
Time 7.41, L2-error 1.53e+00, Max error 3.00e+00
Time 7.61, L2-error 1.53e+00, Max error 3.00e+00
Time 7.81, L2-error 1.53e+00, Max error 3.00e+00
Time 8.01, L2-error 1.53e+00, Max error 3.00e+00
Time 8.21, L2-error 1.53e+00, Max error 3.00e+00
Time 8.41, L2-error 1.53e+00, Max error 3.00e+00
Time 8.61, L2-error 1.53e+00, Max error 3.00e+00
Time 8.81, L2-error 1.53e+00, Max error 3.00e+00
Time 9.01, L2-error 1.53e+00, Max error 3.00e+00
Time 9.21, L2-error 1.53e+00, Max error 3.00e+00
Time 9.41, L2-error 1.53e+00, Max error 3.00e+00
Time 9.61, L2-error 1.53e+00, Max error 3.00e+00
Time 9.81, L2-error 1.53e+00, Max error 3.00e+00
Time 10.00, L2-error 1.53e+00, Max error 3.00e+00

Topic		Replies	Views
Implementing boundary conditions for slanted edges	5	277	September 29, 2022
Strange behavior with boundary conditions and PETSc solvers in dolfin 2019.2	2	465	April 9, 2021
Gmsh: quadrilateral mesh mesh mesh , dolfinx , gmsh	3	74	April 16, 2024
3D vector field: plot 2D slice of 3D mesh dolfinx mesh	5	88	April 12, 2024
From gmsh to dolfin without file generation mesh	5	1300	February 10, 2022

Using Poiseuille example (Test problem 1: Channel flow) with more complex geometries generated with gmsh and additional questions regarding units, implementation of Power Law for viscosity, and best approach for what I'm trying to accomplish

Related Topics