Dear dokken,
I use dolfin-adjoint development version (quay.io/dolfinadjoint/pyadjoint:latest) and was trying to use mumps solver to solve a linear problem (Stokes flow) in a part of the domain. However, I have got the below error.
Exception has occurred: 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 solve linear system using PETSc Krylov solver.
*** Reason: Solution failed to converge in 0 iterations (PETSc reason DIVERGED_PC_FAILED, residual norm ||r|| = 0.000000e+00).
*** Where: This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: f0a90b8fffc958ed612484cd3465e59994b695f9
*** -------------------------------------------------------------------------
File "/workspaces/old-mixed-dimensional/workspace/dolfin-adjoint/Sensitivity_analysis/Channel_flow/Stokes_flow/dolfin-adjoint/dolfin-adjonit_Thermal_Sensitivity.py", line 200, in <module>
dJdsT = JhatT.derivative()
The Stokes flow problem can be solved successfully when I use the default linear solver, but at some point when the mesh is fine enough, it also return an error UMFPACK V5.7.1 (Oct 10, 2014): ERROR: out of memory (whcih should not be the case bacause I was tracking the memory usage and there was still some more free memory left, e.g. 500 MB)
Interestingly, after solving this Stokes flow problem, I solved thermal conduction and convection equation in the same script which also using the solution from the Strokes flow in the weak form, and I called nonlinear solver to solve this thermal problem, which the linear solver is set to mump. The later problem can be solved without arising of any error.
The below is the MWE,
from fenics import *
from fenics_adjoint import *
# classes for sub domains
class Fluid(SubDomain):
def inside(self, x, on_boundary):
return x[1] >= 0.0
class Solid(SubDomain):
def inside(self, x, on_boundary):
return x[1] <= 0.0
# classes for boundaries
class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.1)
class LeftF(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0) and between(x[1], (0.0, 0.1))
class RightF(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 1.0) and between(x[1], (0.0, 0.1))
class Interface(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.0)
class LeftS(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0) and between(x[1], (-0.1, 0.0))
class RightS(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 1.0) and between(x[1], (-0.1, 0.0))
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], -0.1)
# Fluid parameters for the variational problem
muF = Constant(1.002, name="muF")
rhoF = Constant(998, name="rhoF")
cdF = Constant(0.5984, name="cdF") #0.5984
cpF = Constant(4181, name="cpF")
alpF = Constant(cdF/rhoF/cpF, name="alpF")
# Solid parameters for the variational problem
muS = Constant(1, name="muS")
rhoS = Constant(2705, name="rhoS")
cdS = Constant(227, name="cdS") #227
cpS = Constant(900, name="cpS")
alpS = Constant(cdS/rhoS/cpS, name="alpS")
mesh = RectangleMesh(Point(0.0, -0.1), Point(1.0, 0.1), 1000, 100)
# Define markers for Subdomians, Dirichlet and Neuman bcs
subdomain_marker = MeshFunction("size_t", mesh, mesh.geometry().dim(), 0)
Fluid().mark(subdomain_marker, 1)
Solid().mark(subdomain_marker, 2)
boundary_marker = MeshFunction("size_t", mesh, mesh.geometry().dim()-1, 0)
Top().mark(boundary_marker, 1)
LeftF().mark(boundary_marker, 2)
RightF().mark(boundary_marker, 3)
Interface().mark(boundary_marker, 4) # interface
LeftS().mark(boundary_marker, 5)
RightS().mark(boundary_marker, 6)
Bottom().mark(boundary_marker, 7)
# Define mesh perturbation field
S = VectorFunctionSpace(mesh, "CG", 1)
s = Function(S, name="Mesh perturbation field")
ALE.move(mesh, s)
# Solving fluid problem
# Define Function spaces
V = VectorElement("CG", triangle, 2) # Velocity on the fluid domain
P = FiniteElement("CG", triangle, 1) # Pressure on the fluid domain
element = MixedElement([V, P])
W = FunctionSpace(mesh, element)
# Define boundary conditions for fluid problem
U_noslip = Constant((0.0,0.0), name="U_noslip")
pressure_left = Constant(15.0, name="pressure_left")
pressure_right = Constant(0.0, name="pressure_right")
bc1 = DirichletBC(W.sub(0), U_noslip, boundary_marker, 1)
bc2 = DirichletBC(W.sub(0), U_noslip, boundary_marker, 4)
bcs = [bc1,bc2]
# define trial functions and test functions
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
# define the domains for integration
dx = Measure("dx", domain=mesh, subdomain_data=subdomain_marker)
ds = Measure("ds", domain=mesh, subdomain_data=boundary_marker)
# unit normal vector on the surfaces of the fluid mesh
n = FacetNormal(mesh)
# Linear solver
# weak form for Stroke flow
a = muF*inner(grad(u), grad(v))*dx(1) \
- p*div(v)*dx(1) \
+ q*div(u)*dx(1)
l = - pressure_left*inner(n, v)*ds(2)\
- pressure_right*inner(n, v)*ds(3)
# Assemble the LHS and RHS
A = assemble(a, keep_diagonal= True)
L = assemble(l)
# Apply boundary conditions
[bc.apply(A) for bc in bcs]
[bc.apply(L) for bc in bcs]
# solve for the solution
A.ident_zeros()
w = Function(W, name="fluid variables")
solve(A, w.vector(), L, 'mumps')
# Save the solutions for the fluid problem
(u, p) = w.split(True)
# Solve the heat transfer problem
# Define Function spaces
T = FiniteElement("CG", triangle, 3)
WT = FunctionSpace(mesh, T)
# Define the boundary conditions for the heat transfer problem
Ti = Constant(300.0, name="T_in")
Tb = Constant(310.0, name="T_bottom")
bc4 = DirichletBC(WT, Ti, boundary_marker, 2)
bc5 = DirichletBC(WT, Tb, boundary_marker, 7)
bcT = [bc4, bc5]
# Define trial and test functions
t = TrialFunction(WT)
z = TestFunction(WT)
# define functions (unknowns in the weak form)
t = Function(WT, name="temperature")
# weak form for the thermal problem
T = cdS*inner(grad(t), grad(z))*dx(2) \
+ cdF*inner(grad(t), grad(z))*dx(1) \
+ rhoF*cpF*inner(dot(u,grad(t)),z)*dx(1)\
# Determine the Jacobian
JT = derivative(T, t)
# solve the nonlinear problem
problem = NonlinearVariationalProblem(T, t, bcT, J=JT)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm['newton_solver']['absolute_tolerance'] = 2.40E-8
prm['newton_solver']['relative_tolerance'] = 4.90E-15
prm['newton_solver']['maximum_iterations'] = 100
prm['newton_solver']['relaxation_parameter'] = 1.0
prm["newton_solver"]["linear_solver"] = "mumps"
solver.solve()
# Set the objective function
JThermal = assemble(t*dx)
# Calculate the sensitivity
JhatT = ReducedFunctional(JThermal, Control(s))
dJdsT = JhatT.derivative()
I would like to ask if this is a bug in the dolfin-adjoint or is there somethig wrong in my script?
Thank you very much
Phuris