I am running a simulation on a mesh with large number of nodes. I need to capture the results at a subset of the nodes. I am running into trouble when I execute the job in parallel. Kindly take a look at the lines in red to understand the problem. While I am
not getting an error, I am also not able to get the correct results here. I would appreciate your assistance in this regard.
###########################################
# Obtaining the solution for a subset of mesh nodes
##########################################
x,y,z = SpatialCoordinate(W.mesh())
phi=ufl.mathfunctions.Atan2(y,x)
Nsteps=int(T/deltat)
time=np.linspace(0,T,Nsteps+1)
resnodesx,rescoords=axialslice(mesh,Ro,np.pi/2)
# axialslice returns a set of nodes and their coordinates corresponding to a particular criterion
resnodes=resnodesx.astype(np.int32) #PETSc gives an error if the indices are not int32
# Trying to save the solution components at all time instants for a particular subset
uradial=np.zeros((len(resnodes),len(time)+3))
utang=np.zeros((len(resnodes),len(time)+3))
uaxial=np.zeros((len(resnodes),len(time)+3))
uradial[:,:3]=rescoords
utang[:,:3]=rescoords
uaxial[:,:3]=rescoords
for i in range(len(time)):
t = time[i]
collar_disp.assign([0.0,0.0,collar_load(t,T_pulse,f_c,amp=dispamp)])
solver1.solve()
# Update old fields with new quantities
u_nm1.assign(u_n)
u_n.assign(u_np1)
u_r.assign(interpolate(u_n[0]*cos(phi)+u_n[1]*sin(phi),W))
u_t.assign(interpolate(-u_n[0]*sin(phi)+u_n[1]*cos(phi),W))
u_z.assign(interpolate(u_n[2],W))
with u_r.dat.vec_ro as vur: #Obtained this idea from Printing in Parallel page of Firedrake Docs
urad=vur.getValues(resnodes)
print(len(urad)) # length of urad seems to differ for each process
uradial[:,i+3]=urad
with u_t.dat.vec_ro as vut:
utan=vut.getValues(resnodes)
utang[:,i+3]=utan
with u_z.dat.vec_ro as vuz:
uax=vuz.getValues(resnodes)
uaxial[:,i+3]=uax
if i%20 ==0:
outfile.write(u_n,u_r,u_t,u_z, time=t)
gc.collect()
#Consequently the CSV's contain only a subset of the nodal results
np.savetxt(csvradfilename,uradial,fmt='%1.6e',delimiter=',')
np.savetxt(csvtangfilename,utang,fmt='%1.6e',delimiter=',')
np.savetxt(csvaxialfilename,uaxial,fmt='%1.6e',delimiter=',')