Solving a no linear equation of the Spalart-Allmaras model

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!

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.