Hi everyone,
I currently have a fine cylindrical mesh that I built in Gmsh of the form:
R = 1.5;
lc = .15;
Point(1) = {0,0,0,lc};
Point(2) = {R,0,0,lc};
Point(3) = {0,R,0,lc};
Point(4) = {-R,0,0,lc};
Point(5) = {0,-R,0,lc};
Circle(1) = {2,1,3};
Circle(2) = {3,1,4};
Circle(3) = {4,1,5};
Circle(4) = {5,1,2};
Line Loop(5) = {1,2,3,4};
Plane Surface(6) = {5};
Extrude {0,0,-8} {
Surface{6};
}
Field[1] = Distance;
Field[1].NodesList = {1};
Field[2] = Threshold;
Field[2].IField = 1;
Field[2].LcMin = lc / 30;
Field[2].LcMax = lc;
Field[2].DistMin = 0.01;
Field[2].DistMax = 4;
Field[3] = Min;
Field[3].FieldsList = {2};
Background Field = 3;
Mesh.CharacteristicLengthExtendFromBoundary = 0;
Mesh.CharacteristicLengthFromPoints = 0;
Mesh.CharacteristicLengthFromCurvature = 0;
I have then imported this into my code via
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2)
But when I run it, I get this errror:
fenics@aef3e8f96120:~/shared$ python3 target3Dmsh.py
Time 0.0
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 6.335e-01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-03)
Traceback (most recent call last):
File "target3Dmsh.py", line 86, in <module>
solve(F==0, T, [], J = Gain, solver_parameters={"newton_solver":{"linear_solver": "mumps", "relative_tolerance": 1e-3} }, form_compiler_parameters={"cpp_optimize": True, "representation": "uflacs"} )
File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py", line 220, in solve
_solve_varproblem(*args, **kwargs)
File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py", line 266, in _solve_varproblem
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 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.1.0
*** Git changeset: 74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------
I think that this is an issue with the solver because I am using a very fine mesh, but I am unsure on how I can fix this. My full code is shown:
from __future__ import print_function
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
from dolfin import *
from mshr import *
from ufl import as_tensor
from ufl import Index
import math
import mshr
parameters["allow_extrapolation"] = True
parameters["form_compiler"]["cpp_optimize"] = True
set_log_level(20)
# SOLVE THE HEAT TRANSFER PROBLEM FOR A 3D TARGET
#----------------
t_start = 0.0
t_end = 20
nstep = 40
dtime = (t_end - t_start)/ nstep #time step
#-----------------
#Configure for dimensions of system
#----------THERMAL--PROPERTIES--[Element: W]--------
kappa = 0.0172 #conductivity [W/mm K]
c = 0.132 #Specific Heat Capacity [J/gK]
rho = 0.01925 #Density [g/mm^3]
const = kappa /(c * rho)
tau_T = 0.0 #Temperature Gradient Lag
tau_q = 0.0 #Heat Flux Lag
#alpha = 4.3943e+4 #absorption coefficient (mm^-1)
alpha = 1
#---------------------------------------
pi = 3.141592653589793
T_am = 298 #ambient vacuum temperature (K)
T_a4 = T_am**4
epsilon = 0.25 # material emissivity
sigma = 5.67E-14 # W/(mm**2.K**4)
es = epsilon*sigma
Pmax = 500 #peak power of laser in W
w = 0.5 #mm
R = 1.5 #mm
area = pi*R*R #area of target (mm^2)
depth = 8.0 #thickness of target (mm)
volume = area * depth #volume of target
Iden = Pmax / (pi*w*w*volume*rho) #Intensity per unit mass
absorb_depth = 1/alpha
Laser = Expression('2*A*alpha*Iden*exp(-pow((x[0] - 0), 2)/(w*w)-pow((x[1]-0), 2)/(w*w)-pow((x[2]-0), 2)/(w2*w2))',degree=3, A=0.25, alpha=alpha, Iden=Iden,w= 0.56, w2= absorb_depth) #power (w2 localises the z-coordinates)
#A =! emissivity for system in thermal equillibrium
# Create mesh
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2)
Space = FunctionSpace((mesh), 'P', 1) #define finite element function space, defined via basis functions
VectorSpace = VectorFunctionSpace(mesh, 'P', 1)
cells = MeshFunction('size_t',mesh,mesh.topology().dim()) #codimension of 0
facets = MeshFunction('size_t',mesh,mesh.topology().dim()-1) #codimension of 1
da = Measure('ds', domain=mesh, subdomain_data = facets, metadata = {'quadrature_degree': 2}) #area element
dv = Measure('dx', domain=mesh, subdomain_data = cells, metadata = {'quadrature_degree': 2}) #volume element
initial_T = Expression("Tini", Tini=T_am, degree=3) # extrapolate an expression for the temperature before heating
T0 = interpolate(initial_T, Space)
T = Function(Space) # Temperature
V = TestFunction(Space) # Test Function used for FEA
dT = TrialFunction(Space) # Temperature Derivative
q0 = Function(VectorSpace) # heat flux at previous time step
i = Index()
G = as_tensor(T.dx(i), (i)) #gradient of T
G0 = as_tensor(T0.dx(i), (i)) # gradient of T at previous time step
q = as_tensor(dtime/(dtime + tau_q) * (tau_q/dtime*q0[i] - kappa*(1+tau_T/dtime)*G[i] + kappa*tau_T/dtime*G0[i]),(i)) #heat
F = (rho*c/dtime*(T-T0)*V - q[i]*V.dx(i) - rho*Laser*V ) * dv + es*(T**4 - T_a4)*V*da #final form to solve
Gain = derivative(F, T, dT) # Gain will be usedas the Jacobian required to determine the evolution of a dynamic system
file_T = File('target3D/solution.pvd')
for t in np.arange(t_start,t_end,dtime):
print( "Time", t)
solve(F==0, T, [], J = Gain, solver_parameters={"newton_solver":{"linear_solver": "mumps", "relative_tolerance": 1e-3} }, form_compiler_parameters={"cpp_optimize": True, "representation": "uflacs"} )
file_T << (T,t)
q_tmp = project(q, VectorSpace)
q0.assign(q_tmp)
T0.assign(T) #change so that T0 is equal to T for the next time step