Using Conditional function inside a Project function

Linear Algebra
1

Hi everyone,

I am currently trying to solve a heat equation problem where I want the values of some constants within my heat equation to change when the temperature at a node of my mesh exceed some value (shown as Tmelting in my code). Currently, I am attempting to implement this via a conditional function:

assign(kappa, project(conditional(gt(T, Tmelting), kappa_l, kappa_s), Space, solver_type = "mumps", form_compiler_parameters={"cpp_optimize": True, "representation":"uflacs", "quadrature_degree": 2} ))
# kappa_l = value of kappa if T > Tmelt
# kappa_s = value of kappa if T < Tmelt

Where kappa, kappa_l and kappa_s are all tensors. My aim with this code is to have the value of kappa change from kappa_s to kappa_l when T > Tmelt and have this new value fed back into the weak form for the next time step in my problem. However I am currently getting an error code:

Traceback (most recent call last):
  File "target3D_surface2.py", line 193, in <module>
    assign(kappa, project(conditional(gt(T, Tmelting), kappa_l, kappa_s), Space, solver_type = "mumps", form_compiler_parameters={"cpp_optimize": True, "representation":"uflacs", "quadrature_degree": 2} ))
AttributeError: 'float' object has no attribute '_cpp_object'

Orginially I assumed this was from the floats I had previously used for kappa_l and kappa_s instead of the tensors required by the project function, but now I am now sure. Could anybody be able to help with determining the root of the error?

The relevant code is shown below:

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

parameters["allow_extrapolation"] = True
parameters["form_compiler"]["cpp_optimize"] = True
set_log_level(20)

#Configure for dimensions of system
#----------THERMAL--PROPERTIES--(Solid)--------
kappa_s = 138.0       #conductivity [W/m K] 
c_s = 277.17         #Specific Heat Capacity [J/kg K]
rho_s = 10000.0        #Density [kg/m^3]
tau_T = 0.0 #Temperature Gradient Lag
tau_q = 0.0 #Heat Flux Lag
alpha_s = 5.1162e+7 #absorption coefficient (m^-1)
epsilon_s = 0.05  # material emissivity
#---------------------------------------

#----------THERMAL--PROPERTIES--(liquid)--------
Tmelting = 2896.0
kappa_l = 80.0       #conductivity [W/m K] 
c_l = 350.0         #Specific Heat Capacity [J/kg K]
rho_l = 9330.0        #Density [kg/m^3]
alpha_l = 5.1162e+7 #absorption coefficient (m^-1)
epsilon_l = 0.2  # material emissivity
#---------------------------------------

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
Pmax = 200 #peak power of laser in W
w = 0.56e-3  #m
Iden = Pmax / (pi*w*w)    #Intensity 

Laser = Expression('B*A*Iden*exp(-pow((x[0] - 0), 2)/(w*w)-pow((x[1]-0), 2)/(w*w))',degree=2, A=1, Iden=Iden, w= w, B = 1) #power (w2 localises the z-coordinates)
#----------------------------------------------
data = 'mesh_2mm'
mesh = Mesh()
hdf = HDF5File(mesh.mpi_comm(), data+'.h5', 'r')
hdf.read(mesh, '/mesh', False)
cells = MeshFunction('size_t', mesh, 3)
hdf.read(cells, '/cells')
facets = MeshFunction('size_t', mesh, 2)
hdf.read(facets, '/facets')
#------------------------------------------------

Space = FunctionSpace(mesh, 'P', 1)      #define finite element function space, defined via basis functions
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 

#assign initial values for weak form
kappa = kappa_s      #conductivity [W/m K] 
c = c_s    #Specific Heat Capacity [J/kg K]
rho = rho_s     #Density [kg/m^3]
epsilon = epsilon_s # material emissivity
A = epsilon
es = epsilon*sigma

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)) * dv + es*(T**4 - T_a4)*V*da - A*Laser*V*da(2)   #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-3
solver.parameters["newton_solver"]["absolute_tolerance"] = 1E-4
solver.parameters["newton_solver"]["convergence_criterion"] = "residual"
solver.parameters["newton_solver"]["error_on_nonconvergence"] = True
solver.parameters["newton_solver"]["linear_solver"] = "gmres"
solver.parameters["newton_solver"]["lu_solver"]["symmetric"] = False 
solver.parameters["newton_solver"]["maximum_iterations"] = 100
solver.parameters["newton_solver"]["preconditioner"] = "amg"
solver.parameters["newton_solver"]["relaxation_parameter"] = 1.05
solver.parameters["newton_solver"]["report"] = False

solver.parameters["newton_solver"]["krylov_solver"]["nonzero_initial_guess"] = True 
solver.parameters["newton_solver"]["krylov_solver"]["relative_tolerance"] = 1E-5
solver.parameters["newton_solver"]["krylov_solver"]["absolute_tolerance"] = 1E-6
solver.parameters["newton_solver"]["krylov_solver"]["monitor_convergence"] = False
solver.parameters["newton_solver"]["krylov_solver"]["maximum_iterations"] = 100000

file_T = File('/test/solution.pvd')
t = 0
solver.solve()
kappa = as_tensor(kappa,())
kappa_s = as_tensor(kappa_s,())
kappa_l = as_tensor(kappa_l,())
c = as_tensor(c,())
rho = as_tensor(rho,())
epsilon = as_tensor(epsilon,())
assign(kappa, project(conditional(gt(T, Tmelting), kappa_l, kappa_s), Space, solver_type = "mumps", form_compiler_parameters={"cpp_optimize": True, "representation":"uflacs", "quadrature_degree": 2} ))
file_T << (T,t)
assign(q0, project(q, VectorSpace, solver_type= "bicgstab", form_compiler_parameters={"cpp_optimize": True, "representation": "uflacs", "quadrature_degree": 2} ) )
assign(T0, T)
2

Please add all definitions of kappa and function spaces, such that one can reproduce the error. Please make it a minimal example, removing the time loop and only include those variables needed to produce the error.

This will make it easier for other people to help you, and is of greater help for other people struggling with the same problem.

3

Ok, this has been altered accordingly.

4

The problem is that you try to assign a function (the result of the projection), to a float.
The initial kappa has to be a function.
This is illustrated in the minimal example below. Please note the following:

  1. Your proposed example is not minimal, and contains alot of variables not required to reproduce the error. The example below should be used as a guideline for further inquires.
  2. To represent kappa, I am using a DG space, as kappa is discontinuous. If you decide to use a CG space, you will get transfer regions at the switch of the material.
from dolfin import *
from ufl import as_tensor

kappa_s = 138.0       #conductivity [W/m K]
kappa_l = 80.0       #conductivity [W/m K]

mesh = UnitSquareMesh(10,10)

Space = FunctionSpace(mesh, 'P', 1)
T = project(Expression("x[0]+0.1*x[1]", degree=2), Space)

kappa_space = FunctionSpace(mesh, "DG", 0)
kappa = project(kappa_s,kappa_space)
kappa_s = as_tensor(kappa_s,())
kappa_l = as_tensor(kappa_l,())
Tmelting = 0.5
print(kappa.vector().get_local())
kappa.assign(project(conditional(gt(T, Tmelting), kappa_l, kappa_s), kappa_space))
print(kappa.vector().get_local())
5

Ah I see. Thank you very much! I will attempt this now