MPIRUN - writing out text in parallel

6

Dear @pgnikhil,

There are several things in your code that should be removed for it to be considered as a minimal example. You have several unused variables, commented code and the time loop is not necessary to reproduce your issue.

To solve you problem, you need to realize that

returns a vector of numpy arrays (one per process). Thus, to get the actual average, you need to do the following:


from dolfin import *
import numpy as np

# Scaled variables
L = 1
W = 1
mu = 1
delta = W/L
beta = 1.25
lambda_ = beta

dt = 0.01

# Create mesh and define function space
mesh = BoxMesh(Point(0, 0, 0), Point(L, W, W), 2, 2, 2)
Volume = L*W*W
V = VectorFunctionSpace(mesh, 'CG', 2)
TensorDG0 = TensorFunctionSpace(mesh, 'DG', 0)

# Define boundary condition
tol = 1E-14

Str_Rate = 1.0

def clamped_boundary(x, on_boundary):
    return on_boundary and x[0] < tol

def allNodes(x, on_boundary):       return x[0] > -tol
def X0(x, on_boundary):             return near(x[0],0.0,DOLFIN_EPS)
def Y0(x, on_boundary):             return near(x[1],0.0,DOLFIN_EPS)
def Z0(x, on_boundary):             return near(x[2],0.0,DOLFIN_EPS)

t = 0.0

Displacements = Expression(('x[0]*StrRate*time'),degree = 0,  StrRate = Str_Rate,time = t)

bc1 = DirichletBC(V.sub(0), Displacements, allNodes,method='pointwise')
bc2 = DirichletBC(V.sub(1), Constant(0.), Y0)
bc3 = DirichletBC(V.sub(2), Constant(0.), Z0)

bc = [bc1,bc2,bc3]
# Define strain and stress

def epsilon(u):
    return sym(nabla_grad(u))

def sigma(u):
    return lambda_*div(u)*Identity(d) + 2*mu*epsilon(u)

# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension()  # space dimension
v = TestFunction(V)
f = Constant((0, 0, 0))
T = Constant((0, 0, 0))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds

# Compute solution
u = Function(V)

if MPI.comm_world.rank ==0:
    StrStrOutput = open(f"Output_{MPI.comm_world.size}.txt",'w')
    StrStrOutput.write("strainss[11] \t strainss[22] \t \
                stresses[11] \t stresses[22] ")

    
t+=dt
    
Displacements.time = t

solve(a == L, u, bc)
stress = project(sigma(u),TensorDG0)
strain = project(epsilon(u),TensorDG0)
localstrain = strain.vector().get_local()
localstress = stress.vector().get_local()

MPI.comm_world.barrier()
localstrain = np.hstack(mesh.mpi_comm().allgather(localstrain))
localstress = np.hstack(mesh.mpi_comm().allgather(localstress))
MPI.comm_world.barrier()
if MPI.comm_world.rank ==0:
    StrStrOutput.write("\n"+str(np.average(localstrain[0::9])/Volume)+"\t"+str(np.average(localstrain[4::9])/Volume)+"\t"+str(np.average(localstress[0::9])/Volume)+"\t"+str(np.average(localstress[4::9])/Volume))
    StrStrOutput.flush()