Hi Everyone
I need to incorporate to my code a non-linear piecewise function for Kappa which depend on the temperature.
My code is already running with
def kappa(T):
return alpha0*(0.8383+0.00142*T)**(-0.5853)*(2*E0/2/(1+nu0)+3*E0*nu0/((1+nu0)*(1-2*nu0)))
However if I try to use
def kappa(T):
if T<600:
alpha=alpha0*(0.8383+0.00142*T)**(-0.5853)*(2*E0/2/(1+nu0)+3*E0*nu0/((1+nu0)*(1-2*nu0)))
else:
alpha=alpha0
return alpha
I got this error:
Traceback (most recent call last):
File “non_linear_thermoelasticity_fenics.py”, line 120, in
mech_form = inner(sigma(du, dTheta), eps(u_))*dx - inner(f, u_)dx + qdot(n,u_)ds(104)
File “non_linear_thermoelasticity_fenics.py”, line 115, in sigma
return (lmbdatr(eps(v)) - kappa(Theta)Theta)Identity(2) + 2mueps(v)
File “non_linear_thermoelasticity_fenics.py”, line 67, in kappa
if ge(T,600):
File “/home/numericalmining/anaconda3/envs/fenics2018/lib/python3.7/site-packages/ufl/conditional.py”, line 46, in bool
error(“UFL conditions cannot be evaluated as bool in a Python context.”)
File “/home/numericalmining/anaconda3/envs/fenics2018/lib/python3.7/site-packages/ufl/log.py”, line 172, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: UFL conditions cannot be evaluated as bool in a Python context.
Here is my full code and the mesh to run it (Dropbox - Files.zip - Simplify your life):
from dolfin import *
from mshr import *
from fenics import *
import numpy as np
if has_linear_algebra_backend("Epetra"):
parameters["linear_algebra_backend"] = "Epetra"
has_linear_algebra_backend('PETSc')
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"
# 1. Read mesh
mesh = Mesh("dyke_mesh_super_fine.xml")
facets = MeshFunction("size_t", mesh, "dyke_mesh_super_fine_facet_region.xml")
ds = ds(subdomain_data=facets)
#caracteristica de la malla
#W=6000 #Model width m
#H=3000 #Model height m
#Prof=1500 #Mean depth magma chamber m
#w1=10 #Horizontal length of the magma chamber m
#h1=500 #Vertical length of the magma chamber m
#el (0,0) del modelo está en el vertice superior izquierdo
#2. Read rock properties
data=np.loadtxt('Model_catalog.csv',skiprows=1,delimiter=',',usecols = (2,3,4,5,6,7,8) )
Temperatura=data[0,0]
Moduloelasticidad=data[0,1]
Poisson=data[0,2]
Expansiontermica=data[0,3]
Capacidadcalorifica=data[0,4]
Conductividadtermica=data[0,5]
Presionmagma=data[0,6]
#nmodel=len(Temperatura)
#for i in np.arange(nmodel):
tmax=24*5 #horas
num_step=10
dtime=tmax/num_step
Tintru=Temperatura #grados celsius
rho0=2700 #densidad
E0=Moduloelasticidad #modulo de elasticidad en MPa
nu0=Poisson # razon de poisson, adimensional
alpha0=Expansiontermica #coeficiente de expansion termica, 1/°C
cv0=Capacidadcalorifica#conductividad termica, J/kg°C
k0=Conductividadtermica #conductividad térmica, J/horas/m°C
p0=Presionmagma #Magma pressure, MPa
filename="sol_thermo_Tintro="+str(Tintru)+"_rho="+str(rho0)+"_E0="+str(E0)+"_nu0="+str(nu0)+"_alpha0="+str(alpha0)+"_cv0="+str(cv0)+"_k0="+str(k0)+"_p0="+str(p0)+".xdmf"
Tm=Constant(Tintru) #temperatura magma
T0=Expression('15-0.03*x[1]',degree=1) #gradiente termal corteza
DTdike=Expression('Tm-T0',degree=1,Tm=Tm,T0=T0)
E=Constant(E0)
nu = Constant(nu0)
lmbda = Constant(E*nu/((1+nu)*(1-2*nu)))
mu = Constant(E/2/(1+nu))
rho = Constant(rho0) # density
alpha = Constant(alpha0) # 8e-6 thermal expansion coefficient
#kappa = Constant(alpha*(2*mu + 3*lmbda))
cV = Constant(cv0)*rho # specific heat per unit volume at constant strain
k = Constant(k0) # thermal conductivity
rho_g = 0 # 0 si no hay fuerzas de cuerpo, 0.027 si las hay
f = Constant((0, -rho_g))
q=Constant(p0)
def kappa(T):
return alpha0*(0.8383+0.00142*T)**(-0.5853)*(2*E0/2/(1+nu0)+3*E0*nu0/((1+nu0)*(1-2*nu0)))
#4. Define space of solutions and boundary conditions
Vue = VectorElement('CG', mesh.ufl_cell(), 2) # displacement finite element
Vte = FiniteElement('CG', mesh.ufl_cell(), 1) # temperature finite element
V = FunctionSpace(mesh, MixedElement([Vue, Vte]))
#modelo semi
#bc1 = DirichletBC(V.sub(0).sub(1), Constant(0.), facets,102) #inferior
#bc2 = DirichletBC(V.sub(0).sub(0), Constant(0.), facets,101) #laterales
#bc3 = DirichletBC(V.sub(0).sub(0), Constant(0.), facets,106) #laterales
#bc4 = DirichletBC(V.sub(1), DTdike, facets,104) #interior
#bc5 = DirichletBC(V.sub(1), Constant(0.), facets,101) #inferior
#bc6 = DirichletBC(V.sub(1), Constant(0.), facets,103) #superficie
#bcs = [bc1, bc2, bc3, bc4, bc5, bc6]
#modelo entero
bc1 = DirichletBC(V.sub(0), Constant((0.,0.)), facets,102) #inferior
bc2 = DirichletBC(V.sub(0).sub(0), Constant(0.), facets,101) #laterales
bc4 = DirichletBC(V.sub(1), DTdike, facets,104) #interior
bc5 = DirichletBC(V.sub(1), Constant(0.), facets,102) #inferior
bc6 = DirichletBC(V.sub(1), Constant(0.), facets,103) #superficie
bc7 = DirichletBC(V.sub(1), Constant(0.), facets,101) #laterales
bcs = [bc1, bc2, bc4, bc5, bc6, bc7]
n = FacetNormal(mesh)
#5.Espacios de funciones y formulaciones debiles
U_ = TestFunction(V)
(u_, Theta_) = split(U_)
dU = Function(V)
#dU = TrialFunction(V)
(du, dTheta) = split(dU)
Uold = Function(V)
(uold, Thetaold) = split(Uold)
def eps(v):
return sym(grad(v))
def sigma(v, Theta):
return (lmbda*tr(eps(v)) - kappa(Theta)*Theta)*Identity(2) + 2*mu*eps(v)
dt=Constant(dtime)
T0=project(T0,V.sub(1).collapse())
mech_form = inner(sigma(du, dTheta), eps(u_))*dx - inner(f, u_)*dx + q*dot(n,u_)*ds(104)
therm_form = (cV*(dTheta-Thetaold)/dt*Theta_ +
kappa(dTheta+T0)*T0*tr(eps(du-uold))/dt*Theta_ +
dot(k*grad(dTheta), grad(Theta_)))*dx
form = mech_form + therm_form
U = Function(V)
file_results = XDMFFile(filename)
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
t=0
for i in np.arange(num_step):
t=t+dtime
solve(form == 0, dU, bcs, solver_parameters={"newton_solver":{"relative_tolerance": 1e-6}})
Uold.assign(dU)
u, T = dU.split()
u.rename('Displacement','label')
file_results.write(u, t)
T.rename('Temperature variation','label')
file_results.write(T,t)
Vsig = TensorFunctionSpace(mesh, "DG",degree=1)
sig = Function(Vsig, name="Cauchy_Stress")
sig.assign(project(-sigma(u,T), Vsig))
file_results.write(sig, t)
Epssig = TensorFunctionSpace(mesh, "DG",degree=1)
epsilon = Function(Epssig, name="Strain")
epsilon.assign(project(eps(u), Epssig))
file_results.write(epsilon,t)
sx = sig.sub(0)
sy = sig.sub(3)
txy = sig.sub(1)
def sigma1(sigma_x, sigma_y, tau_xy): # magnitude of first principal stress
return((sigma_x+sigma_y)/2 + sqrt(((sigma_x-sigma_y)/2)**2 + tau_xy**2))
def sigma3(sigma_x, sigma_y, tau_xy): # magnitude of first principal stress
return((sigma_x+sigma_y)/2 - sqrt(((sigma_x-sigma_y)/2)**2 + tau_xy**2))
S1 = sigma1(sx,sy,txy)
S3 = sigma3(sx,sy,txy)
Difstress=S1-S3
Promstress=(S1+S3)/2
Vp = FunctionSpace(mesh, "P",1)
sig1=Function(Vp,name="Sigma 1")
sig1.assign(project(S1,Vp))
sig3=Function(Vp,name="Sigma 3")
sig3.assign(project(S3,Vp))
sigdif=Function(Vp,name="Differential stress")
sigdif.assign(project(Difstress,Vp))
sigprom=Function(Vp,name="Mean Stress")
sigprom.assign(project(Promstress,Vp))
file_results.write(sig1,t)
file_results.write(sig3,t)
file_results.write(sigdif,t)
file_results.write(sigprom,t)
print('Progress: '+str(100*t/tmax)+'%')
Some idea how I should solve it?