I send you the whole code because I don’t know how to do it minimal… Let’s assume that the equations are well implemented.
The final part is where I try to solve the linear system.
import dolfinx
import gmsh
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import math
import ufl
from ufl import (FacetNormal, FiniteElement, Identity, Measure, TestFunction, TrialFunction, VectorElement,
as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym, MixedElement)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,
create_vector, create_matrix, set_bc)
import matplotlib.pyplot as plt
from dolfinx.fem import (Constant, Function, FunctionSpace,
assemble_scalar, dirichletbc, form, locate_dofs_topological, set_bc)
from dolfinx.mesh import create_unit_square
from petsc4py.PETSc import ScalarType
from dolfinx import fem
from dolfinx import geometry
from dolfinx.io import gmshio
import dolfinx.cpp as _cpp
import dolfinx.fem as _fem
from dolfinx.nls.petsc import NewtonSolver
print("Tamaño del comunicador:", MPI.COMM_WORLD.size)
domain, cell_tags, ft = gmshio.read_from_msh("Mesh/cylinder.msh", MPI.COMM_WORLD, 0, gdim=2)
ft.name = "Facet markers"
#Encontrar los id segun las coordenadas
tdim = domain.topology.dim
tree = dolfinx.geometry.bb_tree(domain, tdim)
num_entities_local = domain.topology.index_map(tdim).size_local + domain.topology.index_map(tdim).num_ghosts
entities = np.arange(num_entities_local, dtype = np.int32)
mid_tree = dolfinx.geometry.create_midpoint_tree(domain, tdim,entities)
mesh_geom = domain.geometry.x
tol = 1e-8
def closest_point(mesh, point,tol):
points = point
while True:
try:
entity = dolfinx.geometry.compute_closest_entity(tree, mid_tree, mesh, points)
break
except RuntimeError:
print(points)
for j in range(len(mesh.geometry.x)):
if (np.abs(mesh_geom[j][0]-points[0])< tol and np.abs(mesh_geom[j][1]-points[1])< tol):
return points, j
geom_dof = dolfinx.cpp.mesh.entities_to_geometry(mesh._cpp_object,tdim,entity, False)
mesh_nodes = mesh_geom[geom_dof][0]
dis = dolfinx.geometry.compute_distance_gjk(points, mesh_nodes)
index = -100
for i in range(len(mesh_nodes)):
#print(mesh_nodes[i])
if (np.abs(mesh_nodes[i][0]-points[0]+ dis[0]) < tol and np.abs(mesh_nodes[i][1]-points[1]+ dis[1]) < tol):
index = i
if (index == -100):
return point, index
else:
return points[0] - dis[0], points[1] - dis[1] , geom_dof[0][index]
#Calculo del campo vectorial normal a las superficies
def facet_normal_approximation(V, mt, mt_id, tangent=False):
#jit_options: dict = {}, form_compiler_options: dict = {}):
comm = V.mesh.comm
n = ufl.FacetNormal(V.mesh)
nh = _fem.Function(V)
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
ds = ufl.ds(domain=V.mesh, subdomain_data=mt, subdomain_id=mt_id)
if tangent:
if V.mesh.geometry.dim == 1:
raise ValueError("Tangent not defined for 1D problem")
elif V.mesh.geometry.dim == 2:
a = ufl.inner(u, v) * ds
L = ufl.inner(ufl.as_vector([-n[1], n[0]]), v) * ds
else:
def tangential_proj(u, n):
"""
See for instance:
https://link.springer.com/content/pdf/10.1023/A:1022235512626.pdf
"""
return (ufl.Identity(u.ufl_shape[0]) - ufl.outer(n, n)) * u
c = fem.Constant(V.mesh, [1, 1, 1])
a = ufl.inner(u, v) * ds
L = ufl.inner(tangential_proj(c, n), v) * ds
else:
a = (ufl.inner(u, v) * ds)
L = ufl.inner(n, v) * ds
# Find all dofs that are not boundary dofs
imap = V.dofmap.index_map
all_blocks = np.arange(imap.size_local, dtype=np.int32)
top_blocks = _fem.locate_dofs_topological(V, V.mesh.topology.dim - 1, mt.find(mt_id))
deac_blocks = all_blocks[np.isin(all_blocks, top_blocks, invert=True)]
# Note there should be a better way to do this
# Create sparsity pattern only for constraint + bc
bilinear_form = _fem.form(a)
pattern = _fem.create_sparsity_pattern(bilinear_form)
pattern.insert_diagonal(deac_blocks)
pattern.finalize()
u_0 = _fem.Function(V)
u_0.vector.set(0)
bc_deac = _fem.dirichletbc(u_0, deac_blocks)
A = _cpp.la.petsc.create_matrix(comm, pattern)
A.zeroEntries()
# Assemble the matrix with all entries
form_coeffs = dolfinx.cpp.fem.pack_coefficients(bilinear_form._cpp_object)
form_consts = _cpp.fem.pack_constants(bilinear_form._cpp_object)
_cpp.fem.petsc.assemble_matrix(A, bilinear_form._cpp_object,
form_consts, form_coeffs, [bc_deac._cpp_object])
if bilinear_form.function_spaces[0] is bilinear_form.function_spaces[1]:
A.assemblyBegin(PETSc.Mat.AssemblyType.FLUSH) # type: ignore
A.assemblyEnd(PETSc.Mat.AssemblyType.FLUSH) # type: ignore
_cpp.fem.petsc.insert_diagonal(A, bilinear_form.function_spaces[0], [bc_deac._cpp_object], 1.0)
A.assemble()
linear_form = _fem.form(L)
b = _fem.petsc.assemble_vector(linear_form)
_fem.petsc.apply_lifting(b, [bilinear_form], [[bc_deac]])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) # type: ignore
_fem.petsc.set_bc(b, [bc_deac])
# Solve Linear problem
solver = PETSc.KSP().create(V.mesh.comm) # type: ignore
solver.setType("cg")
solver.rtol = 1e-8
solver.setOperators(A)
solver.solve(b, nh.vector)
nh.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD) # type: ignore
return nh
Lx = 40
def upstream(x):
a = np.less(x[0],0)
return np.logical_and(np.greater(x[0]**2 + x[1]**2, 5**2),a)
def downstream(x):
return np.isclose(x[0],Lx)
def walls(x):
a = np.greater(x[0],0)
b = np.less(x[0],Lx)
return np.logical_and(np.greater(x[0]**2 + x[1]**2, 5**2), np.logical_and(a,b))
def obstacle(x):
return np.less(x[0]**2 + x[1]**2, 5**2)
fdim = domain.topology.dim -1
up_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, upstream)
down_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, downstream)
wall_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, walls)
obstacle_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, obstacle)
########################################################################################################
# Definition of Spaces. (dim and polinomial type)
v_cg1 = VectorElement("CG", domain.ufl_cell(), 1)
s_cg1 = FiniteElement("CG", domain.ufl_cell(), 1)
V = dolfinx.fem.FunctionSpace(domain, v_cg1)
Q = dolfinx.fem.FunctionSpace(domain, s_cg1)
Z = MixedElement([v_cg1, s_cg1, s_cg1])
W = dolfinx.fem.FunctionSpace(domain, Z)
# Boundary subdomains and coditions
########################################################################################################
boundaries = [(1,up_facets),(2,down_facets),(3,wall_facets),(4,obstacle_facets)]
facet_indices, facet_markers = [], []
for (marker, face ) in boundaries:
facet_indices.append(face)
facet_markers.append(np.full_like(face, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = dolfinx.mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)
V1, _ = W.sub(0).collapse()
Q1, _ = W.sub(1).collapse()
Q1_v, _ = W.sub(2).collapse()
#############################################################################################################
# Define ds measure
ds = Measure("ds", domain=domain, subdomain_data=facet_tag)
#############################################################################################################
n = FacetNormal(domain)
nh = facet_normal_approximation(V1, facet_tag, 4, tangent=False)
dataname = 'Data/data.txt'
data_np = np.loadtxt(dataname)
data_aux = data_np[:,0:2]
indices = []
for i in range(len(data_aux)):
a,y,j = closest_point(domain,data_aux[i],0.001)
indices.append(int(j))
u = Function(V)
p = Function(Q)
nu_t1 = Function(Q)
d = Function(Q1)
u_aux = Function(V)
# AnsysFluent to FEniCSX
###################################################################################################################
# 0 = x_coordinate 1 = y_coordinate 2 = nu turbulenta 3 = nu effective 4 = x_velocity 5 = y_velocity 6 = density 7 = pressure #
###################################################################################################################
for i in range(len(indices)):
nu_t1.x.array[indices[i]] = data_np[i,2]
u_aux.x.array[2*indices[i]]=data_np[i,4]
u_aux.x.array[2*indices[i]+1]=data_np[i,5]
d.x.array[indices[i]] = data_np[i,8]
p.x.array[indices[i]]=data_np[i,7]
u.interpolate(u_aux)
nu_t = Function(Q1)
nu_t.interpolate(nu_t1)
# Constantes del modelo
Cv1 = Constant(domain, PETSc.ScalarType(7.1))
nu_p = Constant(domain, PETSc.ScalarType(9.4e-5))
# Lista de valores para nu_t
nu_t_list = nu_t.x.array[:] # Aquí colocas tus valores para nu_t
nu = 9.4e-5 # Valor de la viscosidad cinemática molecular
C_v1 = 7.1
# Función para resolver el polinomio para cada nu_t
def calcular_tilde_nu(nu_t, nu, C_v1):
# Coeficientes del polinomio de cuarto grado
coeffs = [1, -nu_t, 0, 0, -nu_t * (C_v1**3)*nu**3]
# Encontrar las raíces del polinomio
roots = np.roots(coeffs)
# Seleccionar la raíz real positiva
tilde_nu = np.real(roots[roots > 0])
# Selecciona la raíz más pequeña real positiva
return tilde_nu[0] if len(tilde_nu) > 0 else 0.0000000001
# Aplicar la función a cada valor de la lista nu_t
tilde_nu_list = [calcular_tilde_nu(nu_t, nu, C_v1) for nu_t in nu_t_list]
nu_tilde = Function(Q1)
nu_tilde.x.array[:] = tilde_nu_list
nu_tilde.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
grad_u = ufl.grad(u)
Omega = 0.5 * (grad_u - ufl.transpose(grad_u))
V_adj = V # para ψ₍ᵤ₎
Q_adj = Q # para ψ₍ₚ₎
Q1_adj = Q # para ψ₍ν̃₎
# =============================================================================
# 4. Parámetros y expresiones auxiliares (modelo SA con expresiones explícitas)
# =============================================================================
rho_val = 1.225
mu_val = 1e-3
sigma_val = 2.0/3.0
c_b1 = 0.1355
c_b2 = 0.3
c_w1 = 0.3
kappa = 0.41
# --- Relación clásica para la viscosidad turbulenta ---
# Para ρ=1, la viscosidad cinemática ν = mu_val.
# Definimos: χ = nu_tilde/nu = nu_tilde/mu_val.
C_v1 = Constant(domain, PETSc.ScalarType(7.1))
chi_expr = nu_tilde/nu
f_v1_expr = chi_expr**3/(chi_expr**3 + C_v1**3)
f_v2_expr = 1 - chi_expr/(1 + chi_expr * f_v1_expr)
# Derivada de f_v1:
df_v1_dchi = 3 * chi_expr**2 * C_v1**3/(chi_expr**3 + C_v1**3)**2
df_v1_dnu_tilde = df_v1_dchi*(1.0/nu)
df_v2_dnu_tilde = -(((1.0/nu)*(1 + chi_expr*f_v1_expr)
- chi_expr*(1.0/nu)*(f_v1_expr + chi_expr*df_v1_dchi))
/(1+chi_expr*f_v1_expr)**2)
# Definir la viscosidad turbulenta efectiva (primal) de forma explícita:
nu_t_expr = nu_tilde * f_v1_expr
dnu_t_dnu_tilde = f_v1_expr + nu_tilde*(df_v1_dchi/nu)
mu_t_prime = dnu_t_dnu_tilde*rho_val
S_expr = ufl.sqrt(2*ufl.inner(Omega, Omega) + 1e-12)
tilde_S_expr = S_expr + nu_tilde/(kappa**2*(d+1e-4)**2)*f_v2_expr
dtildeS_dnu = 1.0/(kappa**2*(d+1e-4)**2)
# --- Modelo SA: Definir f_w y su derivada ---
# Definir: r = nu_tilde/(tilde_S * kappa^2 * (d+ε)^2)
r_expr = nu_tilde/(tilde_S_expr*kappa**2*(d+1e-4)**2)
C_w2 = Constant(domain, PETSc.ScalarType(0.3))
C_w3 = Constant(domain, PETSc.ScalarType(2.0))
g_expr = r_expr + C_w2*(r_expr**6 - r_expr)
f_w_expr = g_expr * ((1 + C_w3**6)/(g_expr**6 + C_w3**6))**(1/6)
# Derivada de r respecto a nu_tilde:
dr_dnu = 1/(tilde_S_expr*kappa**2*(d+1e-4)**2) - nu_tilde/(tilde_S_expr**2*kappa**2*(d+1e-4)**2)*dtildeS_dnu
# Derivada de f_w respecto a g:
phi_expr = ((1+C_w3**6)/(g_expr**6+C_w3**6))**(1/6)
df_w_dg = phi_expr - g_expr*(6*g_expr**5)/(g_expr**6+C_w3**6)*(1/6)*phi_expr
dg_dr = 1 + C_w2*(6*r_expr**5 - 1)
df_w_dnu = df_w_dg * dg_dr * dr_dnu
# --- Parte difusiva de la ecuación de transporte para ν̃ ---
# Se define: K(ν̃) = μ + σρ ν̃ y K'(ν̃)= σρ.
K_expr = nu + nu_tilde
K_prime = 1
mu_t_expr_total = nu_t_expr*rho_val
# --- Acoplamiento de la dependencia de μ_t en la ecuación de momentum ---
# Usamos la relación explícita definida arriba: mu_t_expr = nu_t_expr.
# Su derivada con respecto a nu_tilde es: dnu_t/dnu_tilde = dnu_t_dnu_tilde.
# =============================================================================
# 5. Medidas para el dominio y la frontera
# =============================================================================
dx = ufl.dx
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tag)
ds_inlet = ds(1)
ds_outlet = ds(2)
ds_walls = ds(3)
ds_cyl = ds(4)
# =============================================================================
# 6. Construir el espacio mixto para el sistema adjunto
# =============================================================================
element_u = V_adj.ufl_element() # para ψ_u (vector)
element_p = Q_adj.ufl_element() # para ψ_p (escalar)
element_nu = Q1_adj.ufl_element() # para ψ_nu (escalares)
mixed_elem = MixedElement([element_u, element_p, element_nu])
W_adj = fem.FunctionSpace(domain, mixed_elem)
# Definir funciones de ensayo y de ensayo en el espacio mixto:
Psi = ufl.TrialFunction(W_adj) # (ψ_u, ψ_p, ψ_nu)
Phi = ufl.TestFunction(W_adj) # (v, q, w)
psi_u_m, psi_p_m, psi_nu_m = ufl.split(Psi)
v, q, w = ufl.split(Phi)
# =============================================================================
# 7. Formular el sistema adjunto monolítico (volumétrico + frontera)
# =============================================================================
# 7.1. Ecuación adjunta de Continuidad:
a_cont = ufl.inner(ufl.div(psi_u_m), q)*dx + ufl.inner(psi_u_m, q*ufl.FacetNormal(domain))*ds_outlet
# 7.2. Ecuación adjunta de Momentum:
a_mom = (
-rho_val*ufl.inner(ufl.dot(u, ufl.grad(psi_u_m)), v)*dx
+ ufl.inner(ufl.transpose(ufl.grad(u))*psi_u_m, v)*dx
- ufl.inner(ufl.grad(psi_p_m), v)*dx
- ufl.inner((mu_val+mu_t_expr_total)*(ufl.grad(psi_u_m)+ufl.transpose(ufl.grad(psi_u_m))), ufl.grad(v))*dx
+ ufl.inner(psi_nu_m*ufl.grad(nu_tilde), v)*dx
+ ufl.inner(ufl.div(c_b1 * nu_tilde * psi_nu_m*(Omega/S_expr)), v)*dx
)
a_mom += ufl.inner(
(mu_val+mu_t_expr_total)*(ufl.grad(psi_u_m)+ufl.transpose(ufl.grad(psi_u_m)))*ufl.FacetNormal(domain)
- rho_val*psi_u_m*ufl.dot(u, ufl.FacetNormal(domain)), v)*ds_outlet
# 7.3. Ecuación adjunta de Transporte (para ν̃):
a_nu_diff = (
(1/sigma_val)*ufl.inner(K_expr*ufl.grad(psi_nu_m), ufl.grad(w))*dx
+ (2*c_b2/sigma_val)*ufl.inner(ufl.grad(psi_nu_m), ufl.grad(w))*dx
+ rho_val*ufl.inner(ufl.grad(nu_tilde), ufl.grad(psi_nu_m))*w*dx
)
a_nu_couple = mu_t_prime*ufl.inner((ufl.grad(u)+ufl.transpose(ufl.grad(u))), ufl.grad(psi_u_m))*w*dx
a_nu_vol = a_nu_diff + a_nu_couple
# Términos de frontera para la parte difusiva (según la derivación):
n = ufl.FacetNormal(domain)
BT_diff = (
- (1/sigma_val)*ufl.inner(psi_nu_m, K_prime*ufl.dot(ufl.grad(nu_tilde), n) + K_expr*ufl.dot(ufl.grad(w), n))*ds_outlet
+ (1/sigma_val)*ufl.inner(K_expr*ufl.dot(ufl.grad(psi_nu_m), n), w)*ds_outlet
- (2*c_b2/sigma_val)*ufl.inner(ufl.dot(ufl.grad(nu_tilde), n), w)*ds_outlet
)
a_nu_total = a_nu_vol + BT_diff
# Términos extra en la ecuación de transporte por producción y destrucción:
a_nu_prod_extra = - c_b1*(tilde_S_expr + nu_tilde*dtildeS_dnu)*psi_nu_m*w*dx
a_nu_dest_extra = c_w1*(f_w_expr + nu_tilde*df_w_dnu)*(nu_tilde/d**2)*psi_nu_m*w*dx
a_nu_total = a_nu_total + a_nu_prod_extra + a_nu_dest_extra
# Sumar todos los bloques del sistema adjunto:
a_total = a_cont + a_mom + a_nu_total
L_total = 0 # Suponemos que no hay fuente externa (o se añaden ∂J/∂variable si se conoce)
# =============================================================================
# 8. Imponer condiciones de frontera para el sistema adjunto
# =============================================================================
V1, _ = W.sub(0).collapse() # Para ψ₍ᵤ₎ (vector)
Q1_sub, _ = W.sub(1).collapse() # Para ψ₍ₚ₎ (escalar)
Q2_sub, _ = W.sub(2).collapse() # Para ψ₍ν̃₎ (escalar)
fdim = domain.topology.dim - 1
def upstream(x):
a = np.less(x[0], 0)
return np.logical_and(np.greater(x[0]**2 + x[1]**2, 5**2), a)
def downstream(x):
return np.isclose(x[0], Lx)
def walls(x):
a = np.greater(x[0], 0)
b = np.less(x[0], Lx)
return np.logical_and(np.greater(x[0]**2 + x[1]**2, 5**2),
np.logical_and(a, b))
def obstacle(x):
return np.less(x[0]**2 + x[1]**2, 5**2)
def all_boundaries(x):
return np.logical_or(np.logical_or(upstream(x), downstream(x)),
np.logical_or(walls(x), obstacle(x)))
# Localizar las entidades de frontera:
up_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, upstream)
wall_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, walls)
cyl_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, obstacle)
down_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, downstream)
all_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim,all_boundaries)
# Usar locate_dofs_topological para cada subespacio:
up_f = dolfinx.fem.locate_dofs_topological((W.sub(0), V1), fdim, up_facets)
wall_f = dolfinx.fem.locate_dofs_topological((W.sub(0), V1), fdim, wall_facets)
cyl_f = dolfinx.fem.locate_dofs_topological((W.sub(0), V1), fdim, cyl_facets)
down_f = dolfinx.fem.locate_dofs_topological((W.sub(1), Q1_sub), fdim, down_facets)
all_f = dolfinx.fem.locate_dofs_topological((W.sub(2), Q1_sub), fdim, all_facets)
# Definir funciones de condición mediante interpolación:
def zero_vector(x):
return np.zeros((2, x.shape[1]), dtype=np.float64)
def cylinder_condition(x):
return np.array([[2.0]*x.shape[1],
[0.0]*x.shape[1]], dtype=np.float64)
u_icond = fem.Function(V1, name="u_inlet_bc")
u_wall = fem.Function(V1, name="u_wall_bc")
u_cyl = fem.Function(V1, name="u_cylinder_bc")
u_icond.interpolate(zero_vector)
u_wall.interpolate(zero_vector)
u_cyl.interpolate(cylinder_condition)
bc_u = dirichletbc(u_icond, up_f, W.sub(0))
bc_w = dirichletbc(u_wall, wall_f, W.sub(0))
bc_cyl = dirichletbc(u_cyl, cyl_f, W.sub(0))
# Para ψₚ (componente 1): se impone 0 en el downstream.
def zero_scalar(x):
return np.zeros((1, x.shape[1]), dtype=np.float64)
p_icond = fem.Function(Q1_sub, name="p_outlet_bc")
p_icond.interpolate(zero_scalar)
bc_p = dirichletbc(p_icond, down_f, W.sub(1))
bc_nu = dirichletbc(p_icond, all_f, W.sub(2))
bc = [bc_u, bc_w, bc_cyl, bc_p, bc_nu]
# =============================================================================
Ar = lhs(a_total)
bl = rhs(a_total)
problem = dolfinx.fem.petsc.LinearProblem(Ar, bl, bcs=bc, petsc_options={"ksp_type": "gmres", "pc_type": "ilu","ksp_max_it": 1000 })
Phi = problem.solve()
psi_u_sol, psi_p_sol, psi_nu_sol = Psi.split()
I can pass on anything you need to understand it. Sorry for the length of the code.
Thank you for your support.