# How to convert equation units if mesh units do not match parameter units

I’m solving Navier-Stokes equations on a 3D cylinder geomertry to simulate a vessel in human. However, the unit of the mesh generated from gmsh seems to be millimeter while the unit of parameters are SI( the unit of density is kg/m^3, the unit of viscosity is kg/(m s) ,the velocity inlet is m/s, the pressure outlet is Pa, gravity is m/s^2)

Below is the code contains equations, how can I convert the equations or the boundary conditions, this problem has been bothering me for many days:

``````dt=0.001
num_steps=100
mu=0.0035
rho=1060
nu=mu/rho
V=VectorFunctionSpace(mesh,"P",3) #velocity space
Q=FunctionSpace(mesh,"P",1) #pressure space
bcu_in=DirichletBC(V,Constant((0.182866,-0.03869,-0.071156)),boundary_markers,in_marker)

bcp_out=DirichletBC(Q,Constant(666.61),boundary_markers,out_marker)
bcu_wall=DirichletBC(V,Constant((0.0,0.0,0.0)),boundary_markers,wall_marker)
bcu=[bcu_in,bcu_wall]
bcp=[bcp_out]
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

# t=n (u,p) unknown
#u_ = Function(V)
#p_ = Function(Q)

# t=n-1 known
#u_1 = Function(V)
#p_1 = Function(Q)
# Create functions
u0 = Function(V)
u1 = Function(V)
p1 = Function(Q)

# Define coefficients
k = Constant(dt)
g = Constant((0, 0,-9.81))

# Tentative velocity step
F1 = (1/k)*inner(u - u0, v)*dx + inner(grad(u0)*u0, v)*dx + \
a1 = lhs(F1)
L1 = rhs(F1)

# Pressure update
L2 = -(rho/k)*div(u1)*q*dx

# Velocity update
a3 = inner(u, v)*dx
L3 = inner(u1, v)*dx - k/rho*inner(grad(p1), v)*dx
``````

You choose the units for your mesh (mm), so you can either: 1. Rescale the mesh to be in meter
2. rescale density, viscosity etc to be in compatible units with a mesh defined in `mm`

I try to rescale the unit of density, viscosity, velocity and pressure, however an error message still appeared.
I changed the density to 1060/(10^9) kg/mm^3, viscosity to 0.0035/(10^3)(kg/ m), pressure *10^-3(outlet boundary condition), velocity *10^3 (mm/s) (inlet boundary condition)and g to 9810 mm/s^2. I can’t find the error of the units, could you help me?

``````dt=0.001
num_steps=100
mu=0.0000035
rho=0.000001060
nu=mu/rho
V=VectorFunctionSpace(mesh,"P",3) #velocity space
Q=FunctionSpace(mesh,"P",1) #pressure space
bcu_in=DirichletBC(V,Constant((182.866,-38.69,-71.156)),boundary_markers,in_marker)

bcp_out=DirichletBC(Q,Constant(0.66661),boundary_markers,out_marker)
bcu_wall=DirichletBC(V,Constant((0.0,0.0,0.0)),boundary_markers,wall_marker)
bcu=[bcu_in,bcu_wall]
bcp=[bcp_out]
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)
``````

And the error message:

``````  File "test.py", line 104, in <module>
solve(A1, u1.vector(), b1, "bicgstab", "default")
File "/public/lib/python3.8/site-packages/dolfin/fem/solving.py", line 227, in solve
return dolfin.la.solver.solve(*args)
File "/public/solver.py", line 72, in solve
return cpp.la.solve(A, x, b, method, preconditioner)
RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to solve linear system using PETSc Krylov solver.
*** Reason:  Solution failed to converge in 0 iterations (PETSc reason DIVERGED_NANORINF, residual norm ||r|| = inf).
*** Where:   This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.1.0
*** Git changeset:  2e001bd1aae8e14d758264f77382245e6eed04b0
*** -------------------------------------------------------------------------

``````

This error is likely because of poor choice of scaling parameters inducing large floating point error in the linear solver, or your choice of parameters would lead to a Reynolds number which yields a convection dominated flow yielding an unstable numerical scheme.

As a rule of thumb, the ideal goal would be to get a standard 2-norm measure of the residual to be in the range of 10^{-3} to 10^3. It will take care on your part to ensure your formulation is correct and appropriately scaled.

1 Like

L and U in the table refers to?
The dimensionless varible should be implemented in which form?
Excuse me I just can’t understand the table. Should I add coefficients to the equation to achieve unit conversion

Hi! I changed the import file to .step and rescale the mesh, and the unit is right now. But it still meet errors. Here is the new code, I think the code of ns equation has no errors, I really don’t know where is the problem…

``````import meshio
import matplotlib.pyplot as plt
from dolfin import *
#boundary face
#boundary cell
mesh = Mesh()
##
mvc = MeshValueCollection("size_t", mesh, 2)
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
boundary_markers=MeshFunction("size_t",mesh,mvc)
##
mvc = MeshValueCollection("size_t", mesh, 3)
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
T=1
in_marker=1
out_marker=2
wall_marker=3
domain_marker=4
#deltat=T/num_steps
dt=0.01
num_steps=100
mu=0.0035
rho=1060
nu=mu/rho
V=VectorFunctionSpace(mesh,"P",3) #velocity space
Q=FunctionSpace(mesh,"P",1) #pressure space
bcu_in=DirichletBC(V,Constant((0.037386,-0.002018,0.196464)),boundary_markers,in_marker)

bcp_out=DirichletBC(Q,Constant(666.61),boundary_markers,out_marker)
bcu_wall=DirichletBC(V,Constant((0.0,0.0,0.0)),boundary_markers,wall_marker)
bcu=[bcu_in,bcu_wall]
bcp=[bcp_out]
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

# t=n (u,p) unknown
#u_ = Function(V)
#p_ = Function(Q)

# t=n-1 known
#u_1 = Function(V)
#p_1 = Function(Q)
# Create functions
u0 = Function(V)
u1 = Function(V)
p1 = Function(Q)

# Define coefficients
k = Constant(dt)
f = Constant((0, 0, 0))

# Tentative velocity step
F1 = (1/k)*inner(u - u0, v)*dx +inner(grad(u0)*u0, v)*dx + \
a1 = lhs(F1)
L1 = rhs(F1)

# Pressure update
L2 = -(rho/k)*div(u1)*q*dx

# Velocity update
a3 = inner(u, v)*dx
L3 = inner(u1, v)*dx - (k/rho)*inner(grad(p1), v)*dx

# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

# Use amg preconditioner if available
prec = "amg" if has_krylov_solver_preconditioner("amg") else "default"

# Use nonzero guesses - essential for CG with non-symmetric BC
parameters['krylov_solver']['nonzero_initial_guess'] = True

# Create files for storing solution

# Time-stepping
t = dt
while t < T + DOLFIN_EPS:
print(t)
# Update pressure boundary condition
#p_in.t = t
begin("Computing tentative velocity")
# Compute tentative velocity step
b1 = assemble(L1)
[bc.apply(A1, b1) for bc in bcu]
solve(A1, u1.vector(), b1, "bicgstab", "default")
end()
begin("Computing pressure correction")
# Pressure correction
b2 = assemble(L2)
[bc.apply(A2, b2) for bc in bcp]
[bc.apply(p1.vector()) for bc in bcp]
solve(A2, p1.vector(), b2, "bicgstab", prec)
end()
begin("Computing velocity correction")
# Velocity correction
b3 = assemble(L3)
[bc.apply(A3, b3) for bc in bcu]
solve(A3, u1.vector(), b3, "bicgstab", "default")
end()
# Save to file
ufile << u1
pfile << p1

# Move to next time step
u0.assign(u1)
t += dt
``````

As we do not have the meshes, it is quite hard to reproduce any behavior that you are seeing.
Please either use a simplified mesh (unit square or cube) or post the full error traceback.

The full script is here, the pipe_ns.py is the solve code. pipe_CAD.step is the geometry. the pipe_msh_cad.py is the mesh feneration script. Thank you very much for helping me exam the script!

The error traceback is the same as the information

``````Traceback (most recent call last):
File "tipe_ns.py", line 104, in <module>
solve(A1, u1.vector(), b1, "bicgstab", "default")
File "/public/home/cxr1/.conda/envs/fenicss/lib/python3.8/site-packages/dolfin/fem/solving.py", line 227, in solve
return dolfin.la.solver.solve(*args)
File "/public/home/cxr1/.conda/envs/fenicss/lib/python3.8/site-packages/dolfin/la/solver.py", line 72, in solve
return cpp.la.solve(A, x, b, method, preconditioner)
RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to solve linear system using PETSc Krylov solver.
*** Reason:  Solution failed to converge in 0 iterations (PETSc reason DIVERGED_NANORINF, residual norm ||r|| = inf).
*** Where:   This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.1.0
*** Git changeset:  2e001bd1aae8e14d758264f77382245e6eed04b0
*** -------------------------------------------------------------------------

``````

for the simplified mesh, I have tested 2D mesh yesterday, and the code seems no problem. Is the problem related to the solver?

I think the problem is either in the grid or in the solver, but I can’t be sure where the problem lies. I’m a bit confused now. I have compared the results and found that the calculated result is correct when the density is 1. However, when the density is changed to 1060, non convergence occurs, and another software based on finite volume can still calculate the result

That makes the Reynolds number 1000 times bigger, i.e. you’ll likely be in a laminar case when density is 1 and a turbulent case when density is 1000. It’s no surprise that the latter case is much more difficult.

If you can solve successfully at low density, then one approach is to iteratively solve using an effective density. Start with the solvable small value and use its solution as an initial guess for the next (nonlinear) iteration, ramping up the density until you reach your target value.

1 Like

The parameters contains density(1060kg/m^3), viscosity(0.0035Pa s) are all properties of blood, and it is assumed as laminar in the finite volume software. So I really can’t understand why it can’t be solved by fenics.

The finite volume scheme is an inherently stabilised scheme as a consequence of interior flux approximation choices. Being low order (first order accurate) it may also be highly dissipative which aids in stability.

The scheme you propose using in your code above is not stabilised. I’m not actually sure if the degree 3 velocity and degree 1 pressure approximations even satisfy the inf-sup condition, or if that’s even necessary for your formulation.

It is not fair to compare two radically different discretisation schemes, one of which is not well suited for the problem you’re seeking to solve.

Finite volume schemes may be implemented with FEniCS also; however, this would require deriving the formulation and its implementation which is quite different from what you have above.

This document gives a helpful review of the fundamental distinction (and similarity) between FEM and FVM methods (and FDM).
https://doi.org/10.1002/9781119176817.ecm2010
FEM centres around the weak formulation, while FVM centres around a balance equation (a conservation law). In certain specific low order conditions they are identical, but in general (in higher order models) they are using a completely different mathematics.

The problem I want to solve is a three dimensional blood vessel geometery, so the degree 3 velocity is necessary for me

Thank you for your suggestion! I will learn the relevant code and make an attempt

The degree of a Finite element space and the dimensionality of a mesh does not need to be the same.

1 Like

If the degree of the Finite element space is 2, and the degree of the mesh is 3, how to define the velocity boundary condition?

Exactly in the same way as you already did

``````bcu_in=DirichletBC(V,Constant((0.037386,-0.002018,0.196464)),boundary_markers,in_marker)
``````

Since your inlet condition is a constant, it is possible to represent it in a polynomial space of any order, including degree 2.