I had a code working fine on version 0.5.0 which stopped to work in the new release 0.6.0, and I cannot find the problem. The subdomains are well defined.
The error:
Traceback (most recent call last):
File "/home/marcel/KCMmodel/FourierIso.py", line 238, in <module>
problem = LinearProblem(a, L, bcs=bc, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
File "/home/marcel/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/dolfinx/fem/petsc.py", line 566, in __init__
self._a = _create_form(a, form_compiler_options=form_compiler_options, jit_options=jit_options)
File "/home/marcel/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/dolfinx/fem/forms.py", line 176, in form
return _create_form(form)
File "/home/marcel/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/dolfinx/fem/forms.py", line 171, in _create_form
return _form(form)
File "/home/marcel/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/dolfinx/fem/forms.py", line 138, in _form
assert all([id(d) == id(data[0]) for d in data])
AssertionError
The code:
import numpy as np
import ufl
from dolfinx.fem.petsc import LinearProblem
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import fem, mesh, plot
from dolfinx.io import XDMFFile
from petsc4py.PETSc import ScalarType
import gmsh
import sys
import argparse
#import system parameters
side = 40.0 #xy size
height = 100.0 #z size
dL = 10 #xy general resolution
dH = 5 #z layer resolution
substrate_res = 20 #z substrate resolution
thickness = [0.06] #Thickness of each layer
c=[2.49*1e-12,2.98*1e-12] #Specific heat c in J/(um)3K
k=[66*1e-6,33*1e-6] #isotropic conductivity in J/(um)2K
tbc = [5e-5] #Thermal boundary conductance of each interface (same order)
pen_depth = 0.012
rpump = 5.0 #Pump beam radius
rprobe = 5.0 #Probe beam radius
beamoffset = 0.0 #Offset between paump and probe beam center
bumppar = -5 #Parameter to increase the resolution in the center
progressionpar = 1.2 #Parameter to increase the resolution close to the surface (1.1-1.3 probably works)
#calculate interface z-coordinates
interfaces = np.cumsum(np.array(thickness))
interfaces = np.append(interfaces,height)
#The frequency "omega" comes from the parameter
parser = argparse.ArgumentParser(description='Calculate phase shift.')
parser.add_argument('freq', type=float, nargs=1,
help='Frequency in Hz')
args = parser.parse_args()
omega = 2*np.pi*args.freq[0]
proc = MPI.COMM_WORLD.rank #For MPI control
# create the mesh
if proc == 0:
gmsh.initialize()
gmsh.option.setNumber("General.Terminal",0)
gmsh.model.add("Model")
lc = 10 #point resolution - not really used
gmsh.model.geo.addPoint(0, 0, 0, lc, 1) #layer 1 points
gmsh.model.geo.addPoint(side, 0, 0, lc, 2)
gmsh.model.geo.addPoint(side, side, 0, lc, 3)
gmsh.model.geo.addPoint(0, side, 0, lc, 4)
gmsh.model.geo.addLine(1, 2, 1) #border layer 1
gmsh.model.geo.mesh.setTransfiniteCurve(1, dL, "Bump", bumppar)
gmsh.model.geo.addLine(2, 3, 2)
gmsh.model.geo.mesh.setTransfiniteCurve(2, dL, "Bump", bumppar)
gmsh.model.geo.addLine(3, 4, 3)
gmsh.model.geo.mesh.setTransfiniteCurve(3, dL, "Bump", bumppar)
gmsh.model.geo.addLine(4, 1, 4)
gmsh.model.geo.mesh.setTransfiniteCurve(4, dL, "Bump", bumppar)
gmsh.model.geo.addCurveLoop([-1, -2, -3, -4], 1)
gmsh.model.geo.addPlaneSurface([1], 1) #1st surface
gmsh.model.geo.mesh.setTransfiniteSurface(1)
gmsh.model.geo.synchronize()
leftwall_list = []
rightwall_list = []
frontwall_list = []
backwall_list = []
interfaces_index = []
for count, i_height in enumerate(interfaces):
gmsh.model.geo.addPoint(0, 0, i_height, lc, 5+4*count)
gmsh.model.geo.addPoint(side, 0, i_height, lc, 6+4*count)
gmsh.model.geo.addPoint(side, side, i_height, lc, 7+4*count)
gmsh.model.geo.addPoint(0, side, i_height, lc, 8+4*count)
gmsh.model.geo.addLine(5+4*count, 6+4*count, 5+8*count)
gmsh.model.geo.mesh.setTransfiniteCurve(5+8*count, dL, "Bump", bumppar)
gmsh.model.geo.addLine(6+4*count, 7+4*count, 6+8*count)
gmsh.model.geo.mesh.setTransfiniteCurve(6+8*count, dL, "Bump", bumppar)
gmsh.model.geo.addLine(7+4*count, 8+4*count, 7+8*count)
gmsh.model.geo.mesh.setTransfiniteCurve(7+8*count, dL, "Bump", bumppar)
gmsh.model.geo.addLine(8+4*count, 5+4*count, 8+8*count)
gmsh.model.geo.mesh.setTransfiniteCurve(8+8*count, dL, "Bump", bumppar)
gmsh.model.geo.addCurveLoop([(5+8*count), (6+8*count), (7+8*count), (8+8*count)], 2+5*count) #top surface edge loops
if count == len(interfaces)-1:
gmsh.model.geo.addLine(1+4*count, 5+4*count, 9+8*count)
gmsh.model.geo.mesh.setTransfiniteCurve(9+8*count, substrate_res, "Progression", progressionpar)
gmsh.model.geo.addLine(2+4*count, 6+4*count, 10+8*count)
gmsh.model.geo.mesh.setTransfiniteCurve(10+8*count, substrate_res, "Progression", progressionpar)
gmsh.model.geo.addLine(3+4*count, 7+4*count, 11+8*count)
gmsh.model.geo.mesh.setTransfiniteCurve(11+8*count, substrate_res, "Progression", progressionpar)
gmsh.model.geo.addLine(4+4*count, 8+4*count, 12+8*count)
gmsh.model.geo.mesh.setTransfiniteCurve(12+8*count, substrate_res, "Progression", progressionpar)
else:
gmsh.model.geo.addLine(1+4*count, 5+4*count, 9+8*count)
gmsh.model.geo.mesh.setTransfiniteCurve(9+8*count, dH)
gmsh.model.geo.addLine(2+4*count, 6+4*count, 10+8*count)
gmsh.model.geo.mesh.setTransfiniteCurve(10+8*count, dH)
gmsh.model.geo.addLine(3+4*count, 7+4*count, 11+8*count)
gmsh.model.geo.mesh.setTransfiniteCurve(11+8*count, dH)
gmsh.model.geo.addLine(4+4*count, 8+4*count, 12+8*count)
gmsh.model.geo.mesh.setTransfiniteCurve(12+8*count, dH)
if count == 0:
gmsh.model.geo.addCurveLoop([1, 10, -5, -9], 3)
gmsh.model.geo.addCurveLoop([2, 11, -6, -10], 4)
gmsh.model.geo.addCurveLoop([3, 12, -7, -11], 5)
gmsh.model.geo.addCurveLoop([4, 9, -8, -12], 6)
else:
gmsh.model.geo.addCurveLoop([(8*count-3), (10+8*count), -(5+8*count), -(9+8*count)], 3+5*count)
gmsh.model.geo.addCurveLoop([(8*count-2), (11+8*count), -(6+8*count), -(10+8*count)], 4+5*count)
gmsh.model.geo.addCurveLoop([(8*count-1), (12+8*count), -(7+8*count), -(11+8*count)], 5+5*count)
gmsh.model.geo.addCurveLoop([(8*count), (9+8*count), -(8+8*count), -(12+8*count)], 6+5*count)
gmsh.model.geo.addPlaneSurface([2+5*count], 6+5*count) #top surface is the last
gmsh.model.geo.mesh.setTransfiniteSurface(6+5*count)
gmsh.model.geo.addPlaneSurface([3+5*count], 3+5*count)
gmsh.model.geo.mesh.setTransfiniteSurface(3+5*count)
gmsh.model.geo.addPlaneSurface([4+5*count], 4+5*count)
gmsh.model.geo.mesh.setTransfiniteSurface(4+5*count)
gmsh.model.geo.addPlaneSurface([5+5*count], 5+5*count)
gmsh.model.geo.mesh.setTransfiniteSurface(5+5*count)
gmsh.model.geo.addPlaneSurface([6+5*count], 2+5*count)
gmsh.model.geo.mesh.setTransfiniteSurface(2+5*count)
lastsurface = 6+5*count #keep track of last surface
interfaces_index.append(6+5*count)
frontwall_list.append(3+5*count)
rightwall_list.append(4+5*count)
backwall_list.append(5+5*count)
leftwall_list.append(2+5*count)
#create layer volume
gmsh.model.geo.addSurfaceLoop([1+5*count, 2+5*count, 3+5*count, 4+5*count, 5+5*count, 6+5*count])
gmsh.model.geo.addVolume([1+count], 1+count)
gmsh.model.geo.mesh.setTransfiniteVolume(1+count)
gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(3, [1+count], 1+count)
interfaces_index.pop(-1)
gmsh.model.addPhysicalGroup(2, [lastsurface], 2)
gmsh.model.addPhysicalGroup(2, [1], 1)
gmsh.model.addPhysicalGroup(2, frontwall_list, 3)
gmsh.model.addPhysicalGroup(2, rightwall_list, 4)
gmsh.model.addPhysicalGroup(2, backwall_list, 5)
gmsh.model.addPhysicalGroup(2, leftwall_list, 6)
for n, i in enumerate(interfaces_index):
gmsh.model.addPhysicalGroup(2, [i], 7+n)
gmsh.model.geo.synchronize()
#gmsh.option.setNumber("Mesh.RecombineAll", 1)
gmsh.model.mesh.generate(3)
gmsh.option.setNumber("Mesh.SaveAll", 0)
gmsh.model.geo.synchronize()
gmsh.write("tempmesh.msh")
#if '-nopopup' not in sys.argv:
#gmsh.fltk.run()
gmsh.finalize()
#read the mesh file
from dolfinx.io.gmshio import read_from_msh
model_rank = 0
domain, cell_tags, facet_tags = read_from_msh("tempmesh.msh", MPI.COMM_WORLD, model_rank, gdim=3)
#Create function space
V = fem.FunctionSpace(domain, ("DG", 1)) #Scalar field works for Fourier
x = ufl.SpatialCoordinate(domain) #spatial coordinates
#define the heat flux at the surface and the volumetric heating
f = 2e-6*((ufl.exp(-x[2]/pen_depth))) * x[0] * (x[0]-side) * x[1] * (x[1]-side) * (2 / (np.pi * rpump * rpump) ) * ufl.exp( - 2 * ((x[0]-side/2.0)**2 + (x[1]-side/2.0)**2) / (rpump * rpump ))
#Variatinol problem
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
#Find boundary points
def boundary_one(x):
return np.logical_or(np.isclose(x[0], 0), np.isclose(x[0],side))
def boundary_two(x):
return np.logical_or(np.isclose(x[1], 0), np.isclose(x[1],side))
def boundary_three(x):
return np.isclose(x[2], height)
#define dirichlet boundary conditions
u_D = fem.Function(V, dtype=np.complex128)
u_D.interpolate(lambda x: 0.0 + 0.0j + x[0]*0 + x[1]*0 + x[2]*0)
bc = []
dofs_D = fem.locate_dofs_geometrical(V, boundary_one)
bc.append(fem.dirichletbc(u_D, dofs_D))
dofs_D = fem.locate_dofs_geometrical(V, boundary_two)
bc.append(fem.dirichletbc(u_D, dofs_D))
dofs_D = fem.locate_dofs_geometrical(V, boundary_three)
bc.append(fem.dirichletbc(u_D, dofs_D))
#Create functions for k and c in the space
Q = fem.FunctionSpace(domain, ("DG", 0))
kappa = fem.Function(Q)
heatcap = fem.Function(Q)
for layer in range(0,len(interfaces)):
foundcells = cell_tags.find(layer+1)
kappa.x.array[foundcells] = np.full_like(foundcells, k[layer], dtype=ScalarType)
heatcap.x.array[foundcells] = np.full_like(foundcells, c[layer], dtype=ScalarType)
#tag surface for selective integration
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags)
dS = ufl.Measure("dS", domain=domain, subdomain_data=facet_tags)
mu = 1 #this is the internal penalty parameter, 1e-2 works well so far
n = ufl.FacetNormal(domain) #normal of each surface
#Fourier problem using SIPM, as 10.1016/j.cam.2006.08.029
a = ( kappa*ufl.inner(ufl.grad(u), ufl.grad(v)) + 1j* omega * heatcap * ufl.inner(u,v) ) * ufl.dx \
- ufl.conj(ufl.inner(ufl.avg(kappa*ufl.grad(v)),ufl.jump(u,n))) * ufl.dS \
- ufl.inner(ufl.avg(kappa*ufl.grad(u)),ufl.jump(v,n))*ufl.dS \
+ mu*ufl.inner(ufl.jump(u),ufl.jump(v)) * ufl.dS
#For each interface is necessary to remove general interface conditions and add the interface boundary conductivity
for i, h in enumerate(tbc):
a -= mu*ufl.inner(ufl.jump(u),ufl.jump(v))*dS(7+i)
a += ufl.inner(ufl.jump(u,n)*h,ufl.jump(v,n))*dS(7+i)
a += ufl.conj(ufl.inner(ufl.avg(kappa*ufl.grad(v)),ufl.jump(u,n)))*dS(7+i)
a += ufl.inner(ufl.avg(kappa*ufl.grad(u)),ufl.jump(v,n))*dS(7+i)
L = + ufl.inner(f, v) * ufl.dx
problem = LinearProblem(a, L, bcs=bc, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
#Export results to a file
import dolfinx.io
uh.name = "Temperature"
kappa.name = "k"
heatcap.name = "c"
with dolfinx.io.XDMFFile(domain.comm, "outputfourier.xdmf", "w") as xdmf:
xdmf.write_mesh(domain)
xdmf.write_function(uh)
xdmf.write_function(kappa)
xdmf.write_function(heatcap)