Solving a no linear equation of the Spalart-Allmaras model

dolfinx
1

Hi everyone!
I want to obtain \overline{\nu} from the no linear equation:
\nu_{t} = \tilde{\nu} \cdot \frac{\left(\frac{\tilde{\nu}}{\nu}\right)^3}{\left(\frac{\tilde{\nu}}{\nu}\right)^3 + C_{v1}^3}

my code is:


Cv1 = Constant(domain, PETSc.ScalarType(7.1))
nu = Constant(domain, PETSc.ScalarType(1e-5))
nu_tilde = dolfinx.fem.Function(Q1)
nu_tilde_test = ufl.TestFunction(Q1)


def chi(nu_tilde):
    return nu_tilde/nu

def fv1(nu_tilde):
    return chi(nu_tilde)**3 / (chi(nu_tilde)**3 + Cv1**3)



F = (nu_t - nu_tilde * fv1(nu_tilde)) * nu_tilde_test * ufl.dx

J = ufl.derivative(F, nu_tilde)                        


problem = dolfinx.fem.petsc.NonlinearProblem(F, nu_tilde, J=J)

from dolfinx.nls.petsc import NewtonSolver

solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
solver.report = True

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
n, converged = solver.solve(nu_tilde)

with dolfinx.io.VTKFile(domain.comm, "Plot_mesh/nu_tilde.vtk", "w") as file:
    file.write_mesh(domain)
    file.write_function(nu_tilde)

and obtain

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[6], line 42
     40 opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
     41 ksp.setFromOptions()
---> 42 n, converged = solver.solve(nu_tilde)
     44 with dolfinx.io.VTKFile(domain.comm, "Plot_mesh/nu_tilde.vtk", "w") as file:
     45     file.write_mesh(domain)

File /usr/lib/petsc/lib/python3/dist-packages/dolfinx/nls/petsc.py:46, in NewtonSolver.solve(self, u)
     43 def solve(self, u: fem.Function):
     44     """Solve non-linear problem into function u. Returns the number
     45     of iterations and if the solver converged."""
---> 46     n, converged = super().solve(u.vector)
     47     u.x.scatter_forward()
     48     return n, converged

RuntimeError: Failed to successfully call PETSc function 'KSPSolve'. PETSc error code is: 75, Arguments are incompatible

Any ideas?

Thank you!

2

Please read Read before posting: How do I get my question answered? . The code you posted isn’t reproducible.

I would first of all try with a direct solver. Are you sure that you can use cg to solve with the jacobian matrix on the left-hand side? Looking at the expressions you have in fv1, I wouldn’t vouch that J is symmetric positive definite.

3

sorry, here is my code:

import dolfinx
import gmsh
from mpi4py import MPI
from petsc4py import PETSc
import pyvista
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 geometry
from dolfinx.io import gmshio
import meshio
import dolfinx.cpp as _cpp
import dolfinx.fem as _fem

domain, cell_tags, ft = gmshio.read_from_msh("Mesh/venturi_b_beta08_60.msh", MPI.COMM_WORLD, 0, gdim=2)
ft.name = "Facet markers"


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

    # Definition of Spaces. (dim and polinomial type)

v_cg1 = VectorElement("CG", domain.ufl_cell(), 1)
v_cg2 = VectorElement("CG", domain.ufl_cell(), 2)
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_cg2, s_cg1, s_cg1])
W = dolfinx.fem.FunctionSpace(domain, Z)



# Constantes del modelo
Cv1 = Constant(domain, PETSc.ScalarType(7.1))
nu = Constant(domain, PETSc.ScalarType(1e-5))
# Aquí debes definir o calcular tu nu_t según tu problema específico
                        
# Definir la variable de Spalart-Allmaras
nu_tilde = dolfinx.fem.Function(Q1)
nu_tilde_test = ufl.TestFunction(Q1)

nu_t = Function(Q1)



# Definir la relación no lineal como un problema variacional
def chi(nu_tilde):
    return nu_tilde/nu

def fv1(nu_tilde):
    return chi(nu_tilde)**3 / (chi(nu_tilde)**3 + Cv1**3)



# Definir la ecuación variacional no lineal
F = (nu_t - nu_tilde * fv1(nu_tilde)) * nu_tilde_test * ufl.dx

# Derivar F respecto a nu_tilde para formar el Jacobiano
J = ufl.derivative(F, nu_tilde)                        

# Crear el problema no lineal
problem = dolfinx.fem.petsc.NonlinearProblem(F, nu_tilde, J=J)

from dolfinx.nls.petsc import NewtonSolver

solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
solver.report = True

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
n, converged = solver.solve(nu_tilde)

with dolfinx.io.VTKFile(domain.comm, "Plot_mesh/nu_tilde.vtk", "w") as file:
    file.write_mesh(domain)
    file.write_function(nu_tilde)

The line of code where the function is loaded is missing, as I am not exactly sure how to do it. I am providing it in various formats along with the mesh in the following link.

What do you mean by a direct solver?

Thankyou

4

Replacing

opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"

with something like

opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
5

I have made the changes and I get the following message:

RuntimeError                              Traceback (most recent call last)
Cell In[12], line 49
     47 opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
     48 ksp.setFromOptions()
---> 49 n, converged = solver.solve(nu_tilde)
     51 with dolfinx.io.VTKFile(domain.comm, "Plot_mesh/nu_tilde.vtk", "w") as file:
     52     file.write_mesh(domain)

File /usr/lib/petsc/lib/python3/dist-packages/dolfinx/nls/petsc.py:46, in NewtonSolver.solve(self, u)
     43 def solve(self, u: fem.Function):
     44     """Solve non-linear problem into function u. Returns the number
     45     of iterations and if the solver converged."""
---> 46     n, converged = super().solve(u.vector)
     47     u.x.scatter_forward()
     48     return n, converged

RuntimeError: Newton solver did not converge because maximum number of iterations reached

I have changed F for testing purposes, even for this simple one


F = (nu_t-nu_tilde**2)*nu_tilde_test*ufl.dx

obtaining the same error.