Hi there,
I am trying to establish an initial guess for the Newton Solver of a non linear variational problem. I have already seen that there are several proposed solutions in this forum that make use of interpolate( Expression() ), but that does not work in my case.
Before posting the code: I am solving a flow inside a pipe with the k-epsilon model. Due to the convergence is really bad, I want to use the solution of a lower Reynolds number as initial guess for the velocities and pressure in the whole problem.
My first thought was to impose the previous solution point by point into the velocities and pressure. When it did not work, I interpolated the previous solutions.
The solver is able to make one iteration, then it stops and this error happens:
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 1.951e+03 (tol = 1.000e-06) r (rel) = 1.000e+00 (tol = 1.000e-10)
Traceback (most recent call last):
File "Opt.py", line 33, in
a= funcs.FlowSolver(mesh,boundary_parts,mu, rho, Um,lbase, Ubase, rhobase,lx,ly)
File "/home/ivan/FEniCSProjects/myPrograms/Turbulento/adimensionalizado con punto inicial 1.2/functions.py", line 148, in FlowSolver
solver.solve()
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
****** fenics-support@googlegroups.com
****** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.
****** -------------------------------------------------------------------------
*** Error: Unable to successfully call PETSc function 'MatSetValuesLocal'.
*** Reason: PETSc error code is: 63.
*** Where: This error was encountered inside /build/dolfin-i8W3lr/dolfin-2016.1.0/dolfin/la/PETScMatrix.cpp.
*** Process: 0
*** *** DOLFIN version: 2016.1.0
*** Git changeset: unknown
*** -------------------------------------------------------------------------
I have tried to implement the same solution as https://fenicsproject.org/qa/3793/jacobian-of-nonlinear-problem-in-mixed-formulation but it does not work.
Any thoughts will be much appreciated.
I leave my code below:
V = VectorFunctionSpace(mesh, 'Lagrange',2)
Q = FunctionSpace(mesh,'CG',1)
initialGuessW=MixedFunctionSpace([V,Q])
initialGuessw=Function(initialGuessW)
initialGuessU,initialGuessP= initialGuess(mesh,boundary_parts,lx,ly,mu,Um,Ubase,initialGuessW,initialGuessw)
W=MixedFunctionSpace([V,Q,Q,Q])
(u, p, k, eps) = w.split(True)
#u=interpolate(initialGuessU,V)
#p=interpolate(initialGuessP,Q)
u.vector().array()[:]=numpy.copy(initialGuessU.vector().array()[:])
p.vector().array()[:]=numpy.copy(initialGuessP.vector().array()[:])
mu_t=C_mu*(k**2.)/eps
Re=rhobase*Ubase*lbase/mu
mu_t_adim=mu_t/(rhobase*Ubase*lbase)
Continuity= q*div(u)*dx
Momentum= (2.*(1./Re)*inner(epsilon(u), epsilon(v)) +2.*mu_t_adim*inner(epsilon(u), epsilon(v)) - div(v)*p +rho*inner(grad(u)*u, v) )*dx
TurbulentEnergy= ( -rho*k*inner(u,grad(th)) -2.*mu_t_adim*inner(epsilon(u),epsilon(u))*th + rho*eps*th + inner(grad(k),grad(th))*(1./mu_t_adim)/sigma_k )*dx
TurbulentDissipation= ( -rho*eps*inner(u, grad(psi)) - 2.*mu_t_adim*eps/k*C_1eps*inner(epsilon(u),epsilon(u))*psi + C_2eps*rho*(eps**2.)/k*psi + inner(grad(psi), grad(eps))*mu_t/sigma_e ) *dx
F = Continuity+Momentum+TurbulentEnergy+TurbulentDissipation
dw = TrialFunction(W)
J = derivative(F,w,dw)
problem = NonlinearVariationalProblem(F,w, bcs, J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm['newton_solver'] ['absolute_tolerance'] = 1E-6
prm['newton_solver'] ['relative_tolerance'] = 1E-10
prm['newton_solver'] ['maximum_iterations'] = 200
prm['newton_solver'] ['relaxation_parameter'] = 0.075
solver.solve()
