# Problem with projection method

**URL:** <https://fenicsproject.discourse.group/t/problem-with-projection-method/17035>\
**Category:** General\
**Created:** [February 20, 2025, 1:24pm UTC](https://fenicsproject.discourse.group/t/problem-with-projection-method/17035 "2025-02-20T13:24:53Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![Karlosjusus](https://avatars.discourse-cdn.com/v4/letter/k/cc9497/32.png) [@Karlosjusus](https://fenicsproject.discourse.group/u/Karlosjusus)\
**Post date:** [February 20, 2025, 1:24pm UTC](https://fenicsproject.discourse.group/t/problem-with-projection-method/17035/1 "2025-02-20T13:24:53Z")

</div>

Hello everyone.  
I am trying to project several functions, obtained from external data. When I project some functions, I get that they are 0, but the original functions are not. In fact there is one function, which does not give zero, but I don’t understand why?

Chi\_expr is equal to zero but S\_expr no.

Any idea?  
In this folder is the mesh and the necessary data for the code to compile.

> **[14.34 MB folder on MEGA](https://mega.nz/folder/RJdjQZJS#ed9c-QbAcsBSht34Z6--FA)**
>
> 2 files

my version is 0.7.3

```auto
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 meshio
import dolfinx.cpp as _cpp
import dolfinx.fem as _fem
from dolfinx.nls.petsc import NewtonSolver

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

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)
    #geom_dof = dolfinx.cpp.mesh.compute_incident_entities(mesh,[entity], 2,1)
    #print('geom',geom_dof)
    mesh_nodes = mesh_geom[geom_dof][0]
    #print('mesh_nodes', mesh_nodes)
    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]

########################################################################################################

    # Definition of Spaces. (dim and polinomial type)

v_cg1 = VectorElement("CG", domain.ufl_cell(), 1)
v_cg2 = VectorElement("CG", domain.ufl_cell(), 1)
s_cg1 = FiniteElement("CG", domain.ufl_cell(), 1)
s_cg2 = FiniteElement("CG", domain.ufl_cell(), 1)
V = dolfinx.fem.FunctionSpace(domain, v_cg2)
Q = dolfinx.fem.FunctionSpace(domain, s_cg2)
Z = MixedElement([v_cg1, s_cg1, s_cg1])
W = dolfinx.fem.FunctionSpace(domain, Z)

    
V1, _ = W.sub(0).collapse()
Q1, _ = W.sub(1).collapse()    

dataname = 'data.txt'
vx_np = np.loadtxt(dataname)
vx_aux = vx_np[:,0:2]

b = []
indices = []
ind = int(0)
for i in range(len(vx_aux)):
    a,y,j = closest_point(domain,vx_aux[i],0.0001)
    indices.append(int(j))

print(len(vx_aux))  

data_np = np.loadtxt(dataname)

u = Function(V1)
nu_t1 = Function(Q)
d = Function(Q)

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]+1]=data_np[i,5]
    
    d.x.array[indices[i]] = data_np[i,8]
    
    u_aux.x.array[2*indices[i]]=data_np[i,4]
    
    
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 None

# 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(Q)
nu_tilde.x.array[:] = tilde_nu_list

grad_u = ufl.grad(u)
Omega = 0.5 * (grad_u - ufl.transpose(grad_u))

S_expr = ufl.sqrt(2*ufl.inner(Omega, Omega))

chi_expr = nu_tilde / nu

def project_expression(expr, space):
    v = ufl.TestFunction(space)
    u_trial = ufl.TrialFunction(space)
    a = ufl.inner(u_trial, v) * ufl.dx
    L = ufl.inner(expr, v) * ufl.dx
    problem = fem.petsc.LinearProblem(a, L)
    return problem.solve()

# Proyectamos cada una de las funciones auxiliares en el espacio S_space:
chi_proj = project_expression(chi_expr, Q)
S_proj = project_expression(S_expr, Q)

def write_function_to_file(func, filename,):
    with dolfinx.io.XDMFFile(MPI.COMM_WORLD, filename, "w") as xdmf:
        xdmf.write_mesh(domain)
        xdmf.write_function(func)

write_function_to_file(chi_proj, "Plot/chi.xdmf")
write_function_to_file(S_proj, "Plot/S.xdmf")

```

---

<div class="post-metadata">

**Author:** ![Stein](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/stein/32/8751_2.png) [@Stein](https://fenicsproject.discourse.group/u/Stein)\
**Post date:** [February 20, 2025, 1:43pm UTC](https://fenicsproject.discourse.group/t/problem-with-projection-method/17035/2 "2025-02-20T13:43:29Z")

</div>

Could you please reduce your code to a _minimal_ working example?

---

<div class="post-metadata">

**Author:** ![Karlosjusus](https://avatars.discourse-cdn.com/v4/letter/k/cc9497/32.png) [@Karlosjusus](https://fenicsproject.discourse.group/u/Karlosjusus)\
**Post date:** [February 20, 2025, 1:53pm UTC](https://fenicsproject.discourse.group/t/problem-with-projection-method/17035/3 "2025-02-20T13:53:17Z")

</div>

I have just reduced it!

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [February 20, 2025, 3:42pm UTC](https://fenicsproject.discourse.group/t/problem-with-projection-method/17035/4 "2025-02-20T15:42:24Z")

</div>

Please reduce it to a **minimal** example. The current example contains way more than a single, problematic variable.

---

<div class="post-metadata">

**Author:** ![Karlosjusus](https://avatars.discourse-cdn.com/v4/letter/k/cc9497/32.png) [@Karlosjusus](https://fenicsproject.discourse.group/u/Karlosjusus)\
**Post date:** [February 20, 2025, 4:14pm UTC](https://fenicsproject.discourse.group/t/problem-with-projection-method/17035/5 "2025-02-20T16:14:35Z")

</div>

This it. I think this is the minimum code that works with my data.

```auto
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, TestFunction, TrialFunction, VectorElement,
                 as_vector, ds, dx, inner, grad, nabla_grad)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector, 
                               create_vector, create_matrix, set_bc)

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

from dolfinx import fem
from dolfinx import geometry
from dolfinx.io import gmshio
import dolfinx.cpp as _cpp
import dolfinx.fem as _fem

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

domain, cell_tags, ft = gmshio.read_from_msh("cylinder.msh", MPI.COMM_WORLD, 0, gdim=2)

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)):
        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]

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)

dataname = '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)
nu_t = Function(Q)

for i in range(len(indices)):
    nu_t.x.array[indices[i]] = data_np[i,2] 
    u.x.array[2*indices[i]+1]=data_np[i,5]
    u.x.array[2*indices[i]]=data_np[i,4]

nu_t_list = nu_t.x.array[:]  

nu = 9.4e-5           

C_v1 = 7.1

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 None

tilde_nu_list = [calcular_tilde_nu(nu_t, nu, C_v1) for nu_t in nu_t_list]

nu_tilde = Function(Q)
nu_tilde.x.array[:] = tilde_nu_list

grad_u = ufl.grad(u)
Omega = 0.5 * (grad_u - ufl.transpose(grad_u))
S_expr = ufl.sqrt(2*ufl.inner(Omega, Omega))
chi_expr = nu_tilde / nu

def project_expression(expr, space):
    v = ufl.TestFunction(space)
    u_trial = ufl.TrialFunction(space)
    a = ufl.inner(u_trial, v) * ufl.dx
    L = ufl.inner(expr, v) * ufl.dx
    problem = fem.petsc.LinearProblem(a, L)
    return problem.solve()

chi_proj = project_expression(chi_expr, Q)
S_proj = project_expression(S_expr, Q)

def write_function_to_file(func, filename,):
    with dolfinx.io.XDMFFile(MPI.COMM_WORLD, filename, "w") as xdmf:
        xdmf.write_mesh(domain)
        xdmf.write_function(func)

write_function_to_file(chi_proj, "Plot/chi.xdmf")
write_function_to_file(S_proj, "Plot/S.xdmf")

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [February 20, 2025, 6:23pm UTC](https://fenicsproject.discourse.group/t/problem-with-projection-method/17035/6 "2025-02-20T18:23:14Z")

</div>

> [@Karlosjusus](#):
>
> Chi\_expr is equal to zero

Start by checking if there are any no-zero values in `nu_tilde`, i.e. what are the values of `tilde_nu_list`?

---

<div class="post-metadata">

**Author:** ![Karlosjusus](https://avatars.discourse-cdn.com/v4/letter/k/cc9497/32.png) [@Karlosjusus](https://fenicsproject.discourse.group/u/Karlosjusus)\
**Post date:** [February 20, 2025, 8:23pm UTC](https://fenicsproject.discourse.group/t/problem-with-projection-method/17035/7 "2025-02-20T20:23:22Z")

</div>

This is a plot of nu\_tilde

 ![nu_tilde](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/9/92880c72701b94d999d7ad2010e782d5ea8a4695.png)

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [February 20, 2025, 8:38pm UTC](https://fenicsproject.discourse.group/t/problem-with-projection-method/17035/8 "2025-02-20T20:38:52Z")

</div>

That is not the `nu_tilde_list` in your code. Printing

```python
print( np.min(nu_tilde.x.array), np.max(nu_tilde.x.array))

```

you get `nan`, thus you have broken something within `calcular_tilde_nu`  
which is clear because

> [@Karlosjusus](#):
>
> ```auto
> 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 None
> 
> ```

where the `return`-type `None` doesn’t make sense. Do you mean `0`?

---

<div class="post-metadata">

**Author:** ![Karlosjusus](https://avatars.discourse-cdn.com/v4/letter/k/cc9497/32.png) [@Karlosjusus](https://fenicsproject.discourse.group/u/Karlosjusus)\
**Post date:** [February 20, 2025, 9:42pm UTC](https://fenicsproject.discourse.group/t/problem-with-projection-method/17035/9 "2025-02-20T21:42:24Z")

</div>

Thank you so much @dokken. These change works, I got nan…! But I have another problem with the proyection. I want to calculate the auxiliary functions of spalart-Allmaras model, (f\_{v1}, f\_{v2}, g, r ,f\_{w})

when I calculate g, I get

 ![r](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/d/dd8c53ade77c8caa462c068f945d37edeedef04b.png)

but  
r=nu\_tilde / (tilde\_S\_expr \* kappa **2 \* d** 2)

and the dominator is

 ![nume](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/e/ee6ee9fe8f8f061c6db667d14654aab06620addc.png)

What do you think is the problem?

And I have another question… I want this variable to obtain the production, diffusion and destruction terms of the spalart model. And I want the derivate of these terms respect velocity and nu\_tilde.

I read that

```auto
v = TrialFunction(V)
dF_du = ufl.derivate(F,u,v)

```

is the partial defivarive of F respect u. But I don’t know if that works because of the way I’ve defined the functions with the projections, I lose the dependence of the functions

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [February 20, 2025, 9:47pm UTC](https://fenicsproject.discourse.group/t/problem-with-projection-method/17035/10 "2025-02-20T21:47:02Z")

</div>

> [@Karlosjusus](#):
>
> But I have another problem with the proyection. I want to calculate the auxiliary functions of spalart-Allmaras model, (f\_{v1}, f\_{v2}, g, r ,f\_{w})
> 
> when I calculate g, I get

I don’t know how you calculate `g`, so I cannot help you further. Please consider trying to reduce the complexity of your code. My initial pass through the code took some time as it wasn’t very clear where the error was. Please spend some time inspecting array values (as I showed above) and see if you spot anything obvious.

---

<div class="post-metadata">

**Author:** ![Karlosjusus](https://avatars.discourse-cdn.com/v4/letter/k/cc9497/32.png) [@Karlosjusus](https://fenicsproject.discourse.group/u/Karlosjusus)\
**Post date:** [February 20, 2025, 10:01pm UTC](https://fenicsproject.discourse.group/t/problem-with-projection-method/17035/11 "2025-02-20T22:01:03Z")

</div>

sorry,… I wanted to say r.

I printed

```auto
print( np.min(nu_tilde.x.array), np.max(nu_tilde.x.array))
print( np.min(nume_proj.x.array), np.max(nume_proj.x.array))

```

and I get

-36.75452931542573 66.7615070997315  
-0.0014174498206385264 18.84129517756522

```auto
grad_u = ufl.grad(u)
Omega = 0.5 * (grad_u - ufl.transpose(grad_u))

S_expr = ufl.sqrt(2*ufl.inner(Omega, Omega))

chi_expr = Function(Q)
f_v1_expr = chi_expr**3 / (chi_expr**3 + C_v1**3)

f_v2_expr = 1 - chi_expr / (1 + chi_expr * f_v1_expr)

tilde_S_expr = S_expr + nu_tilde / (kappa **2 * d** 2) * f_v2_expr
r_expr = nu_tilde / (tilde_S_expr * kappa **2 * d** 2)

nume=tilde_S_expr * kappa **2 * d** 2

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)

```

where kappa is constant, d is the distance to the object. I obtain a good function with the proyection method for almost auxiliary functions, but r, g  
and f\_w.

I use this part od the code:

```auto
def project_expression(expr, space):
    v = ufl.TestFunction(space)
    u_trial = ufl.TrialFunction(space)
    a = ufl.inner(u_trial, v) * ufl.dx
    L = ufl.inner(expr, v) * ufl.dx
    
    problem = fem.petsc.LinearProblem(a, L)
    return problem.solve()

# Proyectamos cada una de las funciones auxiliares en el espacio S_space:
#chi_proj = project_expression(chi_expr, Q)
#Q3 = dolfinx.fem.FunctionSpace(domain, ("CG", 3))

f_v1_proj = project_expression(f_v1_expr, Q)
f_v2_proj = project_expression(f_v2_expr, Q)
S_proj = project_expression(S_expr, Q)
tilde_S_proj = project_expression(tilde_S_expr, Q)
nume_proj = project_expression(nume, Q)
r_proj = project_expression(r_expr, Q)
g_proj = project_expression(g_expr, Q)

```

to save the functions. I tried with High dimension space for the projection, but it doesn’t work.

---

<div class="post-metadata">

**Author:** ![Karlosjusus](https://avatars.discourse-cdn.com/v4/letter/k/cc9497/32.png) [@Karlosjusus](https://fenicsproject.discourse.group/u/Karlosjusus)\
**Post date:** [February 21, 2025, 2:28pm UTC](https://fenicsproject.discourse.group/t/problem-with-projection-method/17035/12 "2025-02-21T14:28:59Z")

</div>

It is possible that there are also problems with the mesh. I will repeat the calculations with new meshes. Regarding the methods for calculating partial derivatives, does it work like this? Thanks again! @dokken