Okay,
The calculations are running so far. Unfortunately I did encounter a new problem.
When using the following code
from dolfinx import Mesh, geometry, Function, DirichletBC, Constant, solve, FunctionSpace, VectorFunctionSpace
from dolfinx.io import XDMFFile, VTKFile
from dolfinx.fem import assemble
from ufl import TrialFunction, TestFunction, dot, grad, dx, ds, inner, Measure, nabla_grad
from mpi4py import MPI
import meshio
import pygmsh as pg
lam = [0.195, 3, 26e-3]
pth = 5.6
rho = 1.18
zair = 1e-3
esize =[2e-4,.5e-4,.5e-4]
>#Meshing
geom = pg.opencascade.Geometry()
iso = geom.add_rectangle([0, 0, 0.0], 5e-3, 2.5e-3, char_length = esize[0])
air = geom.add_rectangle([0.0, 0, 5e-3], 5e-3, 2.5e-3, char_length = esize[1])
axis = [0,0,5e-3]
isok = geom.extrude(iso, axis)
axis_air = [0,0,zair]
airk = geom.extrude(air, axis_air)
sens = geom.add_rectangle([2.1e-3, 0 , 4.2e-3],0.8e-3, 0.8e-3, char_length = esize[2])
axis_sens = [0,0,0.8e-3]
sensk = geom.extrude(sens,axis_sens)
geom.add_raw_code('Mesh.RecombineAll = 1;')
geom.add_raw_code('Mesh.Recombine3DAll = 1;')
geom.add_raw_code('Mesh.Algorithm = 2;')
geom.add_raw_code('Mesh.Algorithm3D = 10;')
geom.add_raw_code('Mesh.Smoothing = 10;')
geom.add_raw_code('General.NumThreads = 4;')
geom.add_raw_code('Geometry.ExtrudeReturnLateralEntities = 1;')
geom.add_raw_code('BooleanDifference{ Volume{1}; Delete; }{ Volume{3}; }')
geom.add_raw_code('Coherence;')
geom.add_raw_code('Physical Surface(1) = {41};') #top
geom.add_raw_code('Physical Surface(2) = {38};') #outflow
geom.add_raw_code('Physical Surface(3) = {39};') #outside
geom.add_raw_code('Physical Surface(4) = {40};') #inflow
geom.add_raw_code('Physical Surface(5) = {32};') #bot
geom.add_raw_code('Physical Surface(6) = {33};') #out Mat
geom.add_raw_code('Physical Surface(7) = {35};')#in mat
geom.add_raw_code('Physical Surface(8) = {36};') #outside Mat
geom.add_physical_volume(isok[1], label = "Iso" )
geom.add_physical_volume(airk[1], label="air")
geom.add_physical_volume(sensk[1], label="sensor")
mesh = pg.generate_mesh(geom, verbose=False)
tetra_cells = np.vstack(np.array([cells.data for cells in mesh.cells
if cells.type == "tetra"]))
triangle_cells = np.vstack(np.array([cells.data for cells in mesh.cells
if cells.type == "triangle"]))
for key in mesh.cell_data_dict["gmsh:physical"].keys():
if key == "triangle":
triangle_data = mesh.cell_data_dict["gmsh:physical"][key]
elif key == "tetra":
tetra_data = mesh.cell_data_dict["gmsh:physical"][key]
tetra_mesh = meshio.Mesh(points=mesh.points, cells={"tetra": tetra_cells},
cell_data={"volume":[tetra_data]})
triangle_mesh =meshio.Mesh(points=mesh.points,
cells=[("triangle", triangle_cells)],
cell_data={"boundary":[triangle_data]})
meshio.write("mesh.xdmf", tetra_mesh)
meshio.write("mf.xdmf", triangle_mesh)
>#calculating
ptherm = pth/(1000*1.6e-3*0.8e-3*0.8e-3)
#Mesh lesen
mesh = Mesh
with XDMFFile(MPI.COMM_WORLD,"mesh.xdmf", "r") as infile:
mesh = infile.read_mesh(dolfinx.cpp.mesh.GhostMode.none, name="Grid")
vols = infile.read_meshtags(mesh, name = "Grid")
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
with XDMFFile(MPI.COMM_WORLD, "mf.xdmf", "r") as xdmf_infile:
mt = xdmf_infile.read_meshtags(mesh, name="Grid")
element =("DG", 0)
F = FunctionSpace(mesh, element)
element =("CG", 2)
V = FunctionSpace(mesh, element)
#Definde boundary condition
ind = [1,2,3,4,5,6,7,8]
dom =[]
#Ordnet Knoten den Flächen der RB zu
for element in ind:
dom.append(mt.indices[mt.values == element])
u_D = Function(V)
u_D.vector.set(0)
bc = []
#RB schreiben
for element in dom:
bc.append(dolfinx.DirichletBC(u_D, dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim-1, element)))
body = []
ind =[9,10,11]
#Knoten und deren Koord holen
x = F.tabulate_dof_coordinates()
cond = Function(F)
vel = Function(F)
f = Function(F)
for i in range(x.shape[0]):
if vols.values[i] == 9:
cond.vector.setValueLocal(i,lam[0])
f.vector.setValueLocal(i,0)
elif vols.values[i] == 11:
cond.vector.setValueLocal(i,lam[1])
f.vector.setValueLocal(i,ptherm)
else:
cond.vector.setValueLocal(i,lam[2])
f.vector.setValueLocal(i,0)
#Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
>a = cond*inner(grad(u), grad(v))*dx
L = f*v*dx
#Compute solution
u = Function(V)
solve(a==L, u, bc)
#save
m = u*dx2
t = assemble.assemble_scalar(m)
overtemp = t/(.8e-3*0.8e-3*0.8e-3)
vtkfile = VTKFile('wuek_automata.pvd')
vtkfile.write(u)
Sorry the code is a bit messy, i chopped huge parts away to simplify it.
Okay now the deal, I need a quite fine mesh, at least for esize[1] and esize[2], something like 0.1e-4. Pc can handle it no problem, but with such small element sizes the calcutation somehow “crashes”.
Crashing means, the calculated temperature is something like 5e-4 instead of lets say 1.5K
I tried several things and if I use a constant thermal conductivity for all bodys this problem does not occur. The problem seems to be a combination of small elements and the huge difference in thermal conductivity betwen “sensor” and “air”.
Is there any way to avoid, whatever is happening there?
Regards
Philipp
