Hi everyone,
I am currently trying to solve a highly non-linear problem, specifically a heat transfer problem. The problem calculates an angle dependent quantity (relative to the normal of the mesh) then uses this in the heat equation.
However I am currently experiencing a very strange problem, if I try use UnitCubeMesh(16,16,16) as my mesh over which I solve my problem. I get an error:
Traceback (most recent call last):
File "target-v4-simplified-2.py", line 103, in <module>
problem = NonlinearVariationalProblem(F, T, [], Gain)
File "/usr/local/fkf/Fenics_openmpi2/Python3/lib64/python3.6/site-packages/dolfin/fem/problem.py", line 157, in __init__
F = Form(F, form_compiler_parameters=form_compiler_parameters)
File "/usr/local/fkf/Fenics_openmpi2/Python3/lib64/python3.6/site-packages/dolfin/fem/form.py", line 22, in __init__
sd = form.subdomain_data()
File "/usr/local/fkf/Fenics_openmpi2/Python3/lib/python3.6/site-packages/ufl/form.py", line 208, in subdomain_data
self._analyze_subdomain_data()
File "/usr/local/fkf/Fenics_openmpi2/Python3/lib/python3.6/site-packages/ufl/form.py", line 438, in _analyze_subdomain_data
if data.ufl_id() != sd.ufl_id():
AttributeError: 'builtin_function_or_method' object has no attribute 'ufl_id'
However, if I use a custom mesh, imported from a .h5 file, I get a completely different error:
File "target-v4-simplified.py", line 127, in <module>
problem = NonlinearVariationalProblem(F, T, [], Gain)
File "/usr/local/fkf/Fenics_openmpi2/Python3/lib64/python3.6/site-packages/dolfin/fem/problem.py", line 157, in __init__
F = Form(F, form_compiler_parameters=form_compiler_parameters)
File "/usr/local/fkf/Fenics_openmpi2/Python3/lib64/python3.6/site-packages/dolfin/fem/form.py", line 44, in __init__
mpi_comm=mesh.mpi_comm())
File "/usr/local/fkf/Fenics_openmpi2/Python3/lib64/python3.6/site-packages/dolfin/jit/jit.py", line 82, in mpi_jit
raise RuntimeError(error_msg)
RuntimeError: Compilation failed on root node.
What is the cause of this? My code is given below:
from __future__ import print_function
from fenics import *
import numpy as np
import scipy as scipy
import matplotlib.pyplot as plt
from dolfin import *
import ufl
from ufl import as_tensor
from ufl import Index
import math
from scipy.optimize import curve_fit
from mpi4py import MPI
from scipy.interpolate import interp1d
from scipy import interp
from scipy.interpolate import UnivariateSpline
parameters["allow_extrapolation"] = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["optimize"] = True
set_log_level(50)
hfont = {'fontname':'Liberation Sans'}
#----------------
t_start = 0.0
t_end = 10
nstep = 10
dtime = (t_end - t_start)/ nstep #time step
#---------------------------------------
pi = 3.141592653589793
T_am = 300 #ambient vacuum temperature (K)
T_a4 = T_am**4
sigma = 5.67E-8 # W/(m**2.K**4) S-B Constant
mesh = UnitCubeMesh(16,16, 16)
n1 = 1.0 #refractive index of vacuum
n2 = 2.9421 # refractive index of W at 515nm
# Calculate the angular dependence of the reflectivity of W
i,j,k = ufl.indices(3)
n = FacetNormal(mesh)
zDir = ufl.as_tensor([0,0,-1]) #incoming direction of light (unit vector)
angle = ufl.as_tensor(acos(n[i] * zDir[i]), ())
Space = FunctionSpace(mesh, 'P', 1) #define finite element function space, defined via basis functions
a,del_a = TrialFunction(Space), TestFunction(Space)
lhsAngle = a*del_a*ds
rhsAngle = angle*del_a*ds
a_proj = Function(Space)
AV = assemble(lhsAngle, keep_diagonal=True)
AV.ident_zeros()
LV = assemble(rhsAngle)
solve(AV, a_proj.vector(), LV)
def R(n1,n2,angle_i):
angle_t = asin((n1/n2) * sin(angle_i))
rs = ((n1*cos(angle_i) - n2*cos(angle_t))/(n1*cos(angle_i) + n2*cos(angle_t)))**2
rp = ((n1*cos(angle_t) - n2*cos(angle_i))/(n1*cos(angle_t) + n2*cos(angle_i)))**2
return 1/2 * (rs + rp)
absorb = 1 - R(n1,n2,a)
x = SpatialCoordinate(mesh)
#--------------------------------------------------
Laser = absorb
VectorSpace = VectorFunctionSpace(mesh, 'P', 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
T = Function(Space)
T0 = Function(Space)
T_init = Expression('Tambient', degree=1, Tambient=300.)
T = project(T_init, Space)
assign(T0,T)
#-------------------------------------------
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 = -1.0*G
F = (1.0/dtime*(T-T0)*V - q[i]*V.dx(i)) * dv + sigma*(T**4 - T_a4)*V*da - Laser*V*da(3) #final form to solve
Gain = derivative(F, T, dT) # Gain will be used as the Jacobian required to determine the evolution of a dynamic system
problem = NonlinearVariationalProblem(F, T, [], Gain)
solver = NonlinearVariationalSolver(problem)
solver.parameters["newton_solver"]["relative_tolerance"] = 1E-4
solver.parameters["newton_solver"]["absolute_tolerance"] = 1E-3
solver.parameters["newton_solver"]["convergence_criterion"] = "residual"
solver.parameters["newton_solver"]["error_on_nonconvergence"] = True
solver.parameters["newton_solver"]["linear_solver"] = "cg"
solver.parameters["newton_solver"]["maximum_iterations"] = 10000
solver.parameters["newton_solver"]["preconditioner"] = "hypre_euclid"
solver.parameters["newton_solver"]["report"] = True
solver.parameters["newton_solver"]["krylov_solver"]["nonzero_initial_guess"] = False
solver.parameters["newton_solver"]["krylov_solver"]["relative_tolerance"] = 1E-4
solver.parameters["newton_solver"]["krylov_solver"]["absolute_tolerance"] = 1E-3
solver.parameters["newton_solver"]["krylov_solver"]["monitor_convergence"] = False
solver.parameters["newton_solver"]["krylov_solver"]["maximum_iterations"] = 100000
t = 0
for t in np.arange(t_start + dtime,t_end + dtime,dtime):
print( "Time", t)
solver.solve()
assign(T0, T)
Note: I am also using UFL version:2019.2.0.dev0