I parallel solve my problem.

```
parameters["ghost_mode"] = "shared_vertex"
A4 = assemble(rig)
```

OR

```
parameters["ghost_mode"] = "shared_facet"
A4 = assemble(rig)
```

But i have this error from parallel.

I parallel solve my problem.

```
parameters["ghost_mode"] = "shared_vertex"
A4 = assemble(rig)
```

OR

```
parameters["ghost_mode"] = "shared_facet"
A4 = assemble(rig)
```

But i have this error from parallel.

Again, you need to produce a minimal working example, as the error you are getting could be because of the order you have your commands in your code.

For your particular problem, my best guess is that you havenâ€™t set the ghost mode parameter **before** creating the mesh.

Following is a minimal code example:

```
from dolfin import *
def test_problem(mode):
parameters["ghost_mode"] = mode
mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, "DG", 1)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(jump(u), jump(v))*dS
A = assemble(a)
if __name__ == "__main__":
modes = ["shared_vertex", "shared_facet", "none"]
for mode in modes:
try:
test_problem(mode)
except RuntimeError:
if MPI.comm_world.rank == 0:
print("Code execution failed for Ghost mode {0:s}".format(mode))
```

For follow-up of new posts, please adhere to the guidelines: Read before posting: How do I get my question answered?

1 Like

I cannot provide a simple solver example that shows this problem. My problem is about mixing two streams. I am solving a 3D problem on a mesh of 600,000 elements. I am getting a memory error when using the LU solver. I tried to solve the problem in parallel, but it led to the above problem. The mumps solver says that the residual norm at the first iteration is zero. I will provide the complete code in the hope that I will get an answer.

```
from dolfin import *
import numpy as np
from mpi4py import MPI
import time
comm = MPI.COMM_WORLD
mesh = Mesh()
with XDMFFile("T.xdmf") as infile:
infile.read(mesh)
V = VectorFunctionSpace(mesh, "CG", 1)
dim = mesh.topology().dim()
if comm.rank == 0:
print('Dim:',dim)
mvc = MeshValueCollection("size_t", mesh, dim-1)
with XDMFFile("facet_mesh_T.xdmf") as infile:
infile.read(mvc, "name_to_read")
facet_domains = cpp.mesh.MeshFunctionSizet(mesh, mvc)
u_n_asd = Function(V)
File("solve/Re1/u_real.xml.gz") >> u_n_asd
F_dg = FunctionSpace(mesh, 'DG', 1)
parameters["ghost_mode"] = "shared_facet"
class DirichletBoundary(SubDomain):
def inside(self, x, on_boundary):
return bool(x[1] > 0.00135)
Dbc_values_1 = Expression('1',degree=1)
bc = DirichletBC(F_dg, Dbc_values_1, DirichletBoundary(), "geometric")
t_n = Function(F_dg)
b = Expression("x[1] > 0.0007 ? 1.0 : 0.0", degree=1)
t_n = interpolate(b, F_dg)
if comm.rank == 0:
print('Define')
t = TrialFunction(F_dg)
ti = TestFunction(F_dg)
Time = 0.01
num_times = 1000
dt = Time/num_times
k = Constant(dt)
alpha = Constant(1e2)
theta = Constant(0.25)
D = Constant(2.63e-10)
f = Constant(0)
# Mesh-related functions
n = FacetNormal(mesh)
h = 2*Circumradius(mesh)
h_avg = (h('+') + h('-'))/2
# ( dot(v, n) + |dot(v, n)| )/2.0
un = (dot(u_n_asd, n) + abs(dot(u_n_asd, n)))/2.0
def a(x,z) :
# Bilinear form
a_int = dot(grad(z), D*grad(x) - u_n_asd*x)*dx
a_fac = D*(alpha/h('+'))*dot(jump(z, n), jump(x, n))*dS \
- D*dot(avg(grad(z)), jump(x, n))*dS \
- D*dot(jump(z, n), avg(grad(x)))*dS
a_vel = dot(jump(z), un('+')*x('+') - un('-')*x('-') )*dS \
+ dot(z, un*x)*ds
a = a_int + a_fac + a_vel
return a
a0=a(t_n,ti)
a1=a(t,ti)
A = (1/dt)*inner(t, ti)*dx - (1/dt)*inner(t_n,ti)*dx + theta*a1 + (1-theta)*a0
F = A
rig = lhs(F)
Lig = rhs(F)
gh = Function(F_dg)
if comm.rank == 0:
print('Assemble')
parameters["ghost_mode"] = "shared_vertex"
A4 = assemble(rig)
bc.apply(A4)
t = 0.0
if comm.rank == 0:
print('Run')
for n in range(num_times):
start_time = time.time()
b4 = assemble(Lig)
bc.apply(b4)
solve(A4, gh.vector(), b4, 'umfpack')
MPI.COMM_WORLD.barrier()
max_t = np.abs(np.array(gh.vector())).max()
max_t_norm = mesh.mpi_comm().allreduce(max_t, op=MPI.MAX)
endtime = time.time() - start_time
if comm.rank == 0:
print('max_u= %.6g endtime = %.2g' % (max_t_norm, endtime))
t_n.assign(gh)
F_cg = FunctionSpace(mesh, 'CG', 1)
up = project(gh, V=F_cg)
File("test_pvd/f.pvd") << up
```

If you had read my reply carefully, you should have seen the following:

I.e. You Need to set the shared facet /shared vertex before loading in your mesh.

There is always a way to reduce the complexity of a code that causes an error. If you do not provide simple examples that can reproduce the error, it becomes extremely challenging for anyone to help you, particularly when you cannot reproduce it using built in meshes (as no one can then reproduce and debug your error).

1 Like

Thanks, this solved my problem! But still I am out of memory, how can I use the MUMPS solver?

```
problem = LinearVariationalProblem(lhs(F),rhs(F), gh, [bc])
solver = LinearVariationalSolver(problem)
...
solver.solve()
```

```
solver = LinearVariationalSolver(problem)
solver.parameters["linear_solver"] = "mumps"
```

To verify this, you can add the line

```
set_log_level(LogLevel.DEBUG)
```

before `solver.solve()`

as it will then print the method used to solve the system

1 Like