Trying to disguise sphere with the Navier-Stokes equations

General
1

good day
I have tried to reproduce the results of [1], but I get the following error:

TypeError                                 Traceback (most recent call last)
File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/function.py:492, in Function.interpolate(self, u, cells, cell_map, expr_mesh, nmm_interpolation_data)
    490 try:
    491     # u is a Function or Expression (or pointer to one)
--> 492     _interpolate(u, cells)
    493 except TypeError:
    494     # u is callable

File /usr/lib/python3.10/functools.py:889, in singledispatch.<locals>.wrapper(*args, **kw)
    886     raise TypeError(f'{funcname} requires at least '
    887                     '1 positional argument')
--> 889 return dispatch(args[0].__class__)(*args, **kw)

File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/function.py:455, in Function.interpolate.<locals>._interpolate(u, cells)
    454 """Interpolate a cpp.fem.Function"""
--> 455 self._cpp_object.interpolate(u, cells, nmm_interpolation_data)

TypeError: interpolate(): incompatible function arguments. The following argument types are supported:
    1. interpolate(self, f: ndarray[dtype=float64, writable=False, shape=(*), order='C'], cells: ndarray[dtype=int32, writable=False, shape=(*), order='C']) -> None
    2. interpolate(self, f: ndarray[dtype=float64, writable=False, shape=(*, *), order='C'], cells: ndarray[dtype=int32, writable=False, shape=(*), order='C']) -> None
    3. interpolate(self, u: dolfinx.cpp.fem.Function_float64, cells: ndarray[dtype=int32, writable=False, shape=(*), order='C'], cell_map: ndarray[dtype=int32, writable=False, shape=(*), order='C'], nmm_interpolation_data: tuple[ndarray[dtype=int32, writable=False, shape=(*), order='C'], ndarray[dtype=int32, writable=False, shape=(*), order='C'], ndarray[dtype=float64, writable=False, shape=(*), order='C'], ndarray[dtype=int32, writable=False, shape=(*), order='C']]) -> None
    4. interpolate(self, expr: dolfinx.cpp.fem.Expression_float64, cells: ndarray[dtype=int32, writable=False, order='C'], expr_mesh: dolfinx.cpp.mesh.Mesh_float64, cell_map: ndarray[dtype=int32, writable=False, order='C']) -> None

Invoked with types: dolfinx.cpp.fem.Function_float64, function, ndarray, dolfinx.fem.function.PointOwnershipData

During handling of the above exception, another exception occurred:

RuntimeError                              Traceback (most recent call last)
Cell In[13], line 60
     57     return values
     59 u_in = fem.Function(V)
---> 60 u_in.interpolate(inflow_profile)
     62 u_zero = fem.Function(V)
     63 u_zero.x.array[:] = 0

File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/function.py:497, in Function.interpolate(self, u, cells, cell_map, expr_mesh, nmm_interpolation_data)
    495 assert callable(u)
    496 x = _cpp.fem.interpolation_coords(self._V.element, self._V.mesh.geometry, cells)
--> 497 self._cpp_object.interpolate(np.asarray(u(x), dtype=self.dtype), cells)

RuntimeError: Interpolation data has the wrong shape/size.

I have tried to update the code with the new updates and this is the error I continue to have. I attach the updated code:

import dolfinx
import numpy as np
from mpi4py import MPI
from dolfinx.cpp.mesh import CellType
from dolfinx.io import (XDMFFile)
from dolfinx.fem import (Constant, DirichletBC, Function, FunctionSpace, apply_lifting, assemble_matrix, 
                         assemble_scalar, assemble_vector, create_vector, locate_dofs_topological, set_bc)
from petsc4py import PETSc
from ufl import (FacetNormal, Identity, Measure, TestFunction, TrialFunction,
                 as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym)
print(f"DOLFINx version: {dolfinx.__version__} is installed")

filename="sphere-4.xdmf"
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "sphere-4.xdmf", "r") as xdmf:
       mesh = xdmf.read_mesh(name="Grid")


#### local mesh ####
## mesh = UnitCubeMesh(MPI.COMM_WORLD, 10, 10, 10)
###################

t = 0.00
T = 1. #8                    # Final time
dt = 1./1600.                # Time step size
num_steps = int(T/dt)
k = Constant(mesh, dt)        
mu = Constant(mesh, 0.001)  # Dynamic viscosity
rho = Constant(mesh, 0.5)     # Density

import dolfinx.fem as fem
Q = fem.functionspace(mesh, ("CG", 1))

eps=1e-14

def Inflow(x):
    return np.abs(x[0] - 0.) < eps

def Outflow(x):
    return np.abs(x[0] - 10.0) < eps

def Walls(x):
    result = np.logical_or(np.abs(x[1] + 0.) < eps, np.abs(x[1] - 5.) < eps)
    result2 = np.logical_or(np.abs(x[2] - 0.) < eps, np.abs(x[2] - 5.) < eps)
    return np.logical_or(result, result2)

def Sphere(x):
    result = np.logical_and(x[0] > 2.0 - eps, x[0] < 4.0 + eps)
    result2 = np.logical_and(x[1] > 2.0 - eps, x[1] < 4.0 + eps)
    result3 = np.logical_and(x[2] > 2.0 - eps, x[2] < 4.0 + eps)
    return np.logical_and(np.logical_and(result, result2), result3)

def inflow_profile(x):
    values = np.zeros((3, x.shape[1]), dtype=PETSc.ScalarType)
    values[0, :] = 1.0
    return values

#def inflow_profile(x):
#    values = np.zeros(x.shape)
#    values[0,:]=1
#    return values
V = fem.functionspace(mesh, ("CG", 2))

u_in = fem.Function(V)
u_in.interpolate(inflow_profile)

u_zero = fem.Function(V)
u_zero.x.array[:] = 0
p_zero = fem.Function(Q)
p_zero.x.array[:] = 0

print("Forma de x:", x.shape)
print("Forma de valores:", values.shape)

inflow_boundary_dofs = fem.locate_dofs_geometrical(V, Inflow)
outflow_boundary_dofs = fem.locate_dofs_geometrical(Q, Outflow)
walls_boundary_dofs = fem.locate_dofs_geometrical(V, Walls)
sphere_boundary_dofs = fem.locate_dofs_geometrical(V, Sphere)

bcu_inflow = fem.DirichletBC(u_in, inflow_boundary_dofs)
bcp_outflow = fem.DirichletBC(p_zero, outflow_boundary_dofs)
bcu_walls = fem.DirichletBC(u_zero, walls_boundary_dofs)
bcu_sphere = fem.DirichletBC(u_zero, sphere_boundary_dofs)

bcu = [bcu_inflow, bcu_walls, bcu_sphere]
bcp = [bcp_outflow]

print(f"DOLFINx version: {dolfinx.__version__} is installed")

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

# Define functions for solutions at previous and current time steps
u_n = Function(V)
u_ = Function(V)
p_n = Function(Q)
p_ = Function(Q)

# Define expressions used in variational forms
U = 0.5 * (u_n + u)
n = FacetNormal(mesh)
f = Constant(mesh, PETSc.ScalarType((0,0,0)))
k = Constant(mesh, PETSc.ScalarType(dt))
mu = Constant(mesh, PETSc.ScalarType(.001))
rho = Constant(mesh, PETSc.ScalarType(1))

# Define symmetric gradient
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 variational problem for step 1
F1 = rho * dot((u - u_n) / k, v) * dx \
     + rho * dot(dot(u_n, nabla_grad(u_n)), v) * dx \
     + inner(sigma(U, p_n), epsilon(v)) * dx \
     + dot(p_n * n, v) * ds - dot(mu * nabla_grad(U) * n, v) * ds \
     - dot(f, v) * dx
a1 = lhs(F1)
L1 = rhs(F1)

# Define variational problem for step 2
a2 = dot(nabla_grad(p), nabla_grad(q)) * dx
L2 = dot(nabla_grad(p_n), nabla_grad(q)) * dx - (1 / k) * div(u_) * q * dx

# Define variational problem for step 3
a3 = dot(u, v) * dx
L3 = dot(u_, v) * dx - k * dot(nabla_grad(p_ - p_n), v) * dx

# Assemble matrices
#A1 = assemble(a1)
#A2 = assemble(a2)
#A3 = assemble(a3)
A1 = assemble_matrix(a1, bcs=bcu)
A1.assemble()
b1 = create_vector(L1)


A2 = assemble_matrix(a2, bcs=bcp)
A2.assemble()
b2 = create_vector(L2)


A3 = assemble_matrix(a3)
A3.assemble()
b3 = create_vector(L3)



print(f"DOLFINx version: {dolfinx.__version__} is installed")



# Solver for step 1
solver1 = PETSc.KSP().create(MPI.COMM_WORLD)
solver1.setOperators(A1)
solver1.setType(PETSc.KSP.Type.BCGS)
pc1 = solver1.getPC()
pc1.setType(PETSc.PC.Type.JACOBI)

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

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


print(f"DOLFINx version: {dolfinx.__version__} is installed")

import dolfinx.io
xdmf = dolfinx.io.XDMFFile(MPI.COMM_WORLD, "solution-4.xdmf", "w")
xdmf.write_mesh(mesh)
xdmf.write_function(u_n, t)
xdmf.write_function(p_n, t)

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
    xdmf.write_function(u_n, t)
    xdmf.write_function(p_n, t)

# Close xmdf file
xdmf.close()

Can you help me with this error, please, I have already tried several solutions but none have worked. Also, I attach files from the old version with .h5 and .xdmf files [2]
DOLFINx version: 0.8.0 is installed

2

The issue is in defining the function space, as you are here defining a second order scalar element, while you want a vector element:

V = fem.functionspace(mesh, ("CG", 2, (mesh.geometry.dim, )))

u_in = fem.Function(V)
u_in.interpolate(inflow_profile)

Please also note that you need to change

to

bcu_inflow = fem.dirichletbc(u_in, inflow_boundary_dofs)
bcp_outflow = fem.dirichletbc(p_zero, outflow_boundary_dofs)
bcu_walls = fem.dirichletbc(u_zero, walls_boundary_dofs)
bcu_sphere = fem.dirichletbc(u_zero, sphere_boundary_dofs)

and

to



# Define variational problem for step 1
F1 = rho * dot((u - u_n) / k, v) * dx \
     + rho * dot(dot(u_n, nabla_grad(u_n)), v) * dx \
     + inner(sigma(U, p_n), epsilon(v)) * dx \
     + dot(p_n * n, v) * ds - dot(mu * nabla_grad(U) * n, v) * ds \
     - dot(f, v) * dx
a1 = dolfinx.fem.form(lhs(F1))
L1 = dolfinx.fem.form(rhs(F1))

# Define variational problem for step 2
a2 = dolfinx.fem.form(dot(nabla_grad(p), nabla_grad(q)) * dx)
L2 = dolfinx.fem.form(dot(nabla_grad(p_n), nabla_grad(q)) * dx - (1 / k) * div(u_) * q * dx)

# Define variational problem for step 3
a3 = dolfinx.fem.form(dot(u, v) * dx)
L3 = dolfinx.fem.form(dot(u_, v) * dx - k * dot(nabla_grad(p_ - p_n), v) * dx)

I would also change the imports

to

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

You should also change the outputting to VTXWriter.
Please note that an up to date splitting scheme for DOLFINx is available at: Test problem 2: Flow past a cylinder (DFG 2D-3 benchmark) — FEniCSx tutorial

1 Like
3

Thank you very much for your help. The program works well, but it consumes a lot of computational resources. Therefore, I wanted to know how to do the same code but now with the steady-state Navier-Stokes equations. If you can help me with a guide, please. Regards.

4

Steady-state Navier-Stokes means that you need to remove the time derivative from your problem. However, you should ensure that the splitting scheme makes sense for this.

The code shouldn’t change too much, as you can simply remove some temporal loops.