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!
