This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

storing function evaluation obtained with parallel DG solver

+1 vote

Now that DG is running in parallel, I'm eager to use it for my problems.

However, I'm facing an issue:
After solving, I want to evaluate and store my solution at a certain point in space at the end of every time step: 

vt = Function(U)
ft = Function(V)

print "actual time stepping"
tstep = 1
t = dt
solus = []

while t < T:
    s.t = t

    # first leap
    b1 = assemble(L1)
    A1solve.solve(vt.vector(), b1)
    q10.assign(vt)

    # second leap
    b2 = assemble(L2)
    A2solve.solve(ft.vector(), b2)
    q20.assign(ft)

    if tstep%100==0:
        elapsed = (time.time() - start_time)
        left = (elapsed/tstep)*(N-tstep)
        print "{:4.0f}/{:4.0f}: {:3.2e} ** est {:4.0f} min left **".format(tstep, N, t, left/60.)
        print "source velocity", vt(S)
    save = vt(R)
    save[0] = t
    solus.append(save)

    tstep += 1
    t = tstep*dt # more accurate than t += dt

np.savetxt("received.txt", np.array(solus))

when running this in parallel I'm using the following parameters

parameters['ghost_mode'] = 'shared_facet' # To be able to run DG in parallel
parameters['allow_extrapolation'] = True # To be able to evaluate my solution

this results in an evaluation on every partition/node, but obviously I'm only interested in the solution from the partition in which the point i'm interested in lies in.

In the end I get one file "received.txt", but this appears to hold only the results of one partition, and probably not the one I want.

The question is:
How can I get the evaluation of every partition in a separate file?
Or, even better, only evaluate in the partition I'm interested in?

Here's a full minimal example:

from dolfin import *
import numpy as np

parameters['ghost_mode'] = 'shared_facet'
parameters['allow_extrapolation'] = True

# Set log level
set_log_level(ERROR)

# The difference function
def dif(l):
    return l('-') - l('+')

# Wave parameters
cl = 6320.  # [m/s]
ct = 3130.  # [m/s]
rho = 2700. # [kg/m^3]
C = cl      # Set max speed for numerical fluxes

# set material parameters (Dummy)
c11 = c22 = c33 = Constant(cl**2*rho) #       [Pa]
c44 = c55 = c66 = Constant(ct**2*rho) #       [Pa]
c12 = c13 = c23 = Constant(c11-2*c66) #       [Pa]

# Define mesh
mesh = RectangleMesh(0., 0., 33.e-3, 27.e-3, 33, 27)
n    = FacetNormal(mesh)

# Define functionspaces
deg = 2
U = VectorFunctionSpace(mesh, "DG", deg, dim = 2)
V = VectorFunctionSpace(mesh, "DG", deg, dim = 3)

# Constants for the simple puls
S = (15e-3, 27e-3) # Source coordinates
R = (10.e-3, 0.) # Receiver coordinates
fc = 6.e5
w  = 3.e-3
a1 = -(pi*fc)**2

# time stepping information
dt = 1.e-8
DT = Constant(dt)
t = dt
T = 2.0325e-7
N = T/dt

# Test and trial functions
q1 = TrialFunction(U) # velocity, deformation
q2 = TrialFunction(V) # velocity, deformation
l1 = TestFunction(U)
l2 = TestFunction(V)

# Define expressions
D = "((.5+a*pow((t-td), 2))*exp(a*pow((t-td), 2)))" # ricker wavelet
apod = "exp(-0.5*pow((x[0]-xs)/w, 2))*(x[1] == ys)"
s = Expression(("0.0", apod+"*"+D), w=w, td=2.e-6, a=a1, xs =S[0], ys = S[1], t = 0.0, degree = 3)
zero2 = Expression(("0.0", "0.0"), degree = 1)
zero3 = Expression(("0.0", "0.0", "0.0"), degree = 1)

# Define selection matrices"
print "building matrices"
Ax1 = as_matrix([[-1., 0. ], [0. , 0. ], [0. , -1.]])
Ax2 = as_matrix([[-c11/rho, -c12/rho, 0], [0, 0, -c66/rho]])
Ay1 = as_matrix([[0, 0], [0, -1], [-1., 0]])
Ay2 = as_matrix([[0, 0, -c66/rho], [-c12/rho, -c22/rho, 0]])

# set inital values
q10 = interpolate(zero2, U)
q20 = interpolate(zero3, V)

# Define fluxes on interior and exterior facets
q1hat    = (n[0]('-')*Ax1 + n[1]('-')*Ay1)*avg(q10) + C*0.5*dif(q20)
q1hatbnd = (n[0]*Ax1 + n[1]*Ay1)*(q10) + C*q20
q2hat    = (n[0]('-')*Ax2 + n[1]('-')*Ay2)*avg(q20) + C*0.5*dif(q10)
q2hatbnd = (n[0]*Ax2 + n[1]*Ay2)*(q20) + C*q10

F1 = inner(q1 - q10, l1)*dx \
  - DT*(inner(Ax2*q20, l1.dx(0)) + inner(Ay2*q20, l1.dx(1)))*dx \
  + DT*inner(q2hat, dif(l1))*dS \
  + DT*inner(q2hatbnd, l1)*ds \
  - DT*inner(s/rho, l1)*dx

F2 = inner(q2 - q20, l2)*dx \
  - DT*(inner(Ax1*q10, l2.dx(0)) + inner(Ay1*q10, l2.dx(1)))*dx \
  + DT*inner(q1hat, dif(l2))*dS \
  + DT*inner(q1hatbnd, l2)*ds

a1, L1 = lhs(F1), rhs(F1)
a2, L2 = lhs(F2), rhs(F2)

A1 = assemble(a1)
A2 = assemble(a2)

# compute LU factorization
A1solve = LUSolver("mumps")
A1solve.set_operator(A1)
A1solve.parameters["reuse_factorization"] = True
A2solve = LUSolver("mumps")
A2solve.set_operator(A2)
A2solve.parameters["reuse_factorization"] = True

# Prepare function for solution
vt = Function(U)
ft = Function(V)

# compute LU factorization
tstep = 1
t = dt
solus = []

# Actual time stepping
while t < T:
    s.t = t

    # first leap
    b1 = assemble(L1)
    A1solve.solve(vt.vector(), b1)
    q10.assign(vt)

    # second leap
    b2 = assemble(L2)
    A2solve.solve(ft.vector(), b2)
    q20.assign(ft)

    solus.append(vt(R))

    tstep += 1
    t = tstep*dt # more accurate than t += dt

np.savetxt("received.txt", np.array(solus))

Thanks in advance!

asked Nov 7, 2014 by stevenvdk FEniCS User (1,590 points)
edited Nov 7, 2014 by stevenvdk

Is this really a problem with DG in parallel? It seems more like a problem of Function
evaluation in parallel. Would it be possible to write a simple Minimal Working Example to illustrate your problem?

commented Nov 7, 2014 by chris_richardson FEniCS Expert (31,740 points)

True,
a general problem when evaluating a solution that was obtained by parallel solve...

I've added a (hopefully) minimal example in the question above.

Is there e.g. a variable I can call to get the id of the current node?

Thanks!

commented Nov 7, 2014 by stevenvdk FEniCS User (1,590 points)

You can get the process number with: MPI.rank(mesh.mpi_comm()) but I'm not sure that's what you need here.
This is more what I would call a MWE:

parameters['allow_extrapolation'] = True
parameters['ghost_mode'] = 'shared_facet' 
mesh = UnitSquareMesh(10,10)

Q = FunctionSpace(mesh, "DG", 1)
F = Function(Q)
F.interpolate(Expression("x[0]*x[1]"))

print F(0.25, 0.25)

When run in parallel, it gives different answers on each process, sometimes way off.
What you need, is a way to determine which process your point is on.

Something like this will work, though not too elegant. If anyone else has a good idea...

import numpy as np

parameters['ghost_mode'] = 'shared_facet'
mesh = UnitSquareMesh(10,10)

Q = FunctionSpace(mesh, "DG", 1)
F = Function(Q)
F.interpolate(Expression("x[0]*x[1]"))

x = np.array([0.25, 0.25])
pt = Point(x)

# Max value of unsigned int                                                                                                            
max_int = 2**32 - 1

rank = MPI.rank(mesh.mpi_comm())
bb = BoundingBoxTree()
bb.build(mesh)
if (bb.compute_first_entity_collision(pt) != max_int):
    print rank, F(x)

Bear in mind, this will print on multiple processes if the cell is ghosted...

commented Nov 7, 2014 by chris_richardson FEniCS Expert (31,740 points)
edited Nov 7, 2014 by chris_richardson

Not too pretty indeed, but it works. Thanks.
Thanks for pointing out my example wasn't minimal, I forgot tho think outside my box ;)

commented Nov 7, 2014 by stevenvdk FEniCS User (1,590 points)

1 Answer

+2 votes

This is one way to get the desired effect. Maybe there should be an issue on bitbucket?

import numpy as np
mesh = UnitSquareMesh(10,10)

Q = FunctionSpace(mesh, "DG", 1)
F = Function(Q)
F.interpolate(Expression("x[0]*x[1]"))

x = np.array([0.25, 0.25])
pt = Point(x)

# Max value of unsigned int                                                                                                             
max_int = 2**32 - 1

rank = MPI.rank(mesh.mpi_comm())
bb = BoundingBoxTree()
bb.build(mesh)
if (bb.compute_first_entity_collision(pt) != max_int):
    print rank, F(x)
answered Nov 7, 2014 by chris_richardson FEniCS Expert (31,740 points)
...