Hi Dokken,
Thanks a lot for your helpful tips. They did help me a lot for which Iâm grateful.
So as you suggested, I separated I/O part and save them in XDMF format. In respect to your previous question about having the XDMF file open somewhere else, I believe that the error was caused by running the XDMFWrite part in parallel because then only one process could do that at a time and the file would be read-only to the other cores. So to solve the problem, I had to split the program into two parts and run the first part in serial just to save the inputs. Then, using mpirun -n 4 python3, I ran the parallel code.
I managed to go through the whole parallel program without any error, but the only problem is that the final output isnât what I expect to get because Iâve already run exactly the same model in serial and I got reasonable answers out of that, but this one isnât even close. I doubt I missed something in the parallel part or did something wrongly. Can you please take a look at the following to snippets in serial and parallel to see if there is a blunder or something causing the problem?
Serial Part:
from fenics import *
import scipy.io as sio
import numpy as np
#%% Constants
Kapa_value = [0.5, 0.55] #Thermal Conductivity (W/m/'c), [Muscle, Tumor] (see âcombine_rawâ m-file)
Cp_value = [3421, 3500]#Heat Capacity (J/Kg/'c), [Muscle, Tumor]
Cb = 3617 #specific heat (Heat Capacity) of blood (J/Kg/'c)
Wb_value = [0.8, 0.5] #local perfusion rate (kg/m3/s), [Muscle, Tumor]
Tb = 37 #Arterial blood temperature
Rho_value = [1090, 1040]#Density of Muscle (Kg/m3)
PwrNormFac = 400; #makeshift (needs to be defined carefully in Matlab)
#%% Modern Mesh conversion using meshio package and xdmf-h5 formats!
import meshio
msh = meshio.read(âGeoGmsh_cylinder.mshâ)
meshio.write(âGeoCylinder.xdmfâ, meshio.Mesh(points=msh.points, cells={âtetraâ: msh.cells[âtetraâ]}))
meshio.write(âcf.xdmfâ, meshio.Mesh(points=msh.points, cells={âtetraâ: msh.cells[âtetraâ]},
cell_data={âtetraâ: {âname_to_readâ:msh.cell_data[âtetraâ][âgmsh:physicalâ]}}))#Physical
meshio.write(âmf.xdmfâ, meshio.Mesh(points=msh.points, cells={âtriangleâ: msh.cells[âtriangleâ]},
cell_data={âtriangleâ: {âname_to_readâ: msh.cell_data[âtriangleâ][âgmsh:physicalâ]}})) #Facet
mesh = Mesh()
with XDMFFile(âGeoCylinder.xdmfâ) as infile:
infile.read(mesh)
mvc = MeshValueCollection(âsize_tâ,mesh,3)
with XDMFFile(mesh.mpi_comm(),âcf.xdmfâ) as infile:
infile.read(mvc, âname_to_readâ)
Markers = cpp.mesh.MeshFunctionSizet(mesh,mvc)
#%% Impressed Source Definition from PLD on FEM nodes!
PLDonFEM = sio.loadmat(âPLD_On_FEM_Nodes_cylinder_650MHz.matâ)
PLD = PLDonFEM[âPLDExportListâ][0]/PwrNormFac
FTermSpace = FunctionSpace(mesh, âCGâ, 1) # so that dofs are only in mesh vertices
f = Function(FTermSpace)
f.vector()[:] = PLD[dof_to_vertex_map(FTermSpace)]
#Write in XDMF format so that it can be retrieved in parallel
PLDXDMF = XDMFFile(âQ.xdmfâ)
PLDXDMF.write_checkpoint(f,âSourceâ,0) #0:for the 1st frequency (in a multi-frequency paradigm>> 1: 2nd freq, âŚ)
#%% Subdomainsâ property assignment
V0 = FunctionSpace(mesh, âDGâ, 0)
Kapa = Function(V0)
Cp = Function(V0)
Wb = Function(V0)
Rho = Function(V0)
for n, cell_no in enumerate(Markers.array()):
Subdomain_Tag = int(cell_no - 1)
Kapa.vector()[n] = Kapa_value[Subdomain_Tag]
Cp.vector()[n] = Cp_value[Subdomain_Tag]
Wb.vector()[n] = Wb_value[Subdomain_Tag]
Rho.vector()[n] = Rho_value[Subdomain_Tag]
#Save them for visualization purposes
MeshExport = File(âCylinder/mesh.pvdâ)
MeshExport << mesh
ThermalCond = File(âCylinder/Kapa.pvdâ)
ThermalCond << Kapa
#Save them for parallelization purposes
KapaXDMF = XDMFFile(âKapa.xdmfâ)
KapaXDMF.write_checkpoint(Kapa,âKâ,0)
CpXDMF = XDMFFile(âCp.xdmfâ)
CpXDMF.write_checkpoint(Cp,âCâ,0)
WbXDMF = XDMFFile(âWb.xdmfâ)
WbXDMF.write_checkpoint(Wb,âWâ,0)
RhoXDMF = XDMFFile(âRho.xdmfâ)
RhoXDMF.write_checkpoint(Rho,âRâ,0)
Parallel Part:
from fenics import *
import scipy.io as sio
import numpy as np
#%%
Cb = 3617 #specific heat (Heat Capacity) of blood (J/Kg/'c)
Tb = 37 #Arterial blood temperature
Bolus_Temp = Constant(20)
Body_Temp = Constant(37) #Normal body temperature which is used as the initial condition throughout the body.
#%% Time Settings
TFrame = 3600 #60 min
Num_Steps = 60
dt = TFrame / Num_Steps #For every second
#%% Modern Mesh conversion
mesh = Mesh()
with XDMFFile(mesh.mpi_comm(),âGeoCylinder.xdmfâ) as infile:
infile.read(mesh)
mvc = MeshValueCollection(âsize_tâ,mesh,3)
with XDMFFile(mesh.mpi_comm(),âcf.xdmfâ) as infile:
infile.read(mvc, âname_to_readâ)
Markers = cpp.mesh.MeshFunctionSizet(mesh,mvc)
mvc = MeshValueCollection(âsize_tâ,mesh,2)
with XDMFFile(mesh.mpi_comm(),âmf.xdmfâ) as infile:
infile.read(mvc, âname_to_readâ)
BCs = cpp.mesh.MeshFunctionSizet(mesh,mvc)
dx = Measure(âdxâ, domain=mesh, subdomain_data=Markers)
#%% Boundary Condition and initial value
V = FunctionSpace(mesh, âCGâ, 2)
Outer_bc = DirichletBC(V, Bolus_Temp, BCs, 3) #DomainBoundary())
T_n = interpolate(Body_Temp, V) #or project(Body_Temp, V)
#%% Define Variational Problem
T = TrialFunction(V)
w = TestFunction(V)
#%% Impressed Source Definition from PLD on FEM nodes!
FTermSpace = FunctionSpace(mesh, âCGâ, 1) # so that dofs are only in mesh vertices
f = Function(FTermSpace)
#read the process-specific portion of the source back from XDMF file stored by the serial part.
with XDMFFile(mesh.mpi_comm(),âQ.xdmfâ) as file:
file.read_checkpoint(f,âSourceâ,0)
#%% Subdomainsâ property assignment
V0 = FunctionSpace(mesh, âDGâ, 0)
Kapa = Function(V0)
Cp = Function(V0)
Wb = Function(V0)
Rho = Function(V0)
#read the subdomain XDMF file stored by the serial part for this process.
with XDMFFile(mesh.mpi_comm(),âKapa.xdmfâ) as file:
file.read_checkpoint(Kapa,âKâ,0)
with XDMFFile(mesh.mpi_comm(),âCp.xdmfâ) as file:
file.read_checkpoint(Cp,âCâ,0)
with XDMFFile(mesh.mpi_comm(),âWb.xdmfâ) as file:
file.read_checkpoint(Wb,âWâ,0)
with XDMFFile(mesh.mpi_comm(),âRho.xdmfâ) as file:
file.read_checkpoint(Rho,âRâ,0)
#%% Symbolic Bioheat equation
F = Twdx + (Kapadt/(CpRho))dot(grad(T), grad(w))dx - (T_n + (dt/(CpRho))f - (CbWbdt/(CpRho))(T - Tb))wdx
a, L = lhs(F), rhs(F)
Create VTK file for saving solution
vtkfile = File(âCylinder/Temperature.pvdâ)
Export back to Matlab
hdfHandle = HDF5File(mesh.mpi_comm(), âCylinder/FinalTemperature.h5â,âwâ)
hdfHandle.write(mesh, â/meshâ)
Time-stepping
T = Function(V)
t = 0
for n in range(Num_Steps):
# Update current time
t += dt
# Compute solution
solve(a == L, T, Outer_bc)
# Save to file and plot solution
vtkfile << (T, t)
# Save in h5 format
Name = '/Temperature/%d_min'%n
hdfHandle.write(T, Name)
# Update previous solution
T_n.assign(T)
BTW, I still have some doubts about where to put mesh.mpi_comm() and where to suppress it in the parallel part. Can the problem possibly come from there?
Thanks a lot again!
