Projection function takes too long, how to calculate components separately?

Linear Algebra
1

Summary: My fenics.project() function takes too long to run and kills the program when my mesh is large enough. I think it can be made more efficient but I am not sure how.

I am using fenics to solve an electrostatics problem and solving for potential. Next, I want to extract the electric field vectors from the potential. My understanding is the best way to do this is

fenics.project(-fenics.grad(PotentialSolutionData))

However this step takes an extremely long time, and the program often fails with “Killed.” I’ve tried using the mumps and other solvers with no improvement (it just tries longer before it fails). I’ve read a similar about a similar problem here, and the user improved the projection function by calculating the projection along each component. However when I tried implementing this user’s solution, I encountered errors. My MWE code is at the end of this file, simply a Unit Cube with wall boundary conditions attempting to solve the Laplacian. The example I found was clearly in a different dimensional space, so I tried adapting to my 3D example, but I keep getting an error that the shapes do not match.

How can I more efficiently calculate electric field from potential data? Would the method of calculating the projection along each component be more efficient in this case? If so, how would I implement it?

I’ve tried a few different dimension, directions, etc. to adapt the previous example, but so far none work, and this is my best guess at what it should be.

import fenics as fn
import numpy

# these are defined previously in code
# sourceSpace = fn.VectorFunctionSpace(mesh, 'CG', 1)
# sourceScalarSpace = fn.FunctionSpace(mesh, 'CG', 1)
# destinationSpace = fn.Function(sourceSpace)
# sourceData = -fn.grad(PotentialSolutionData)


def myProjectFn(sourceSpace, sourceScalarSpace, sourceData, destinationSpace):
    directions = [
        fn.Constant(numpy.array([[1, 0, 0], [0, 0, 0], [0, 0, 0]])) ,
        fn.Constant(numpy.array([[0, 0, 0], [0, 1, 0], [0, 0, 0]])) ,
        fn.Constant(numpy.array([[0, 0, 0], [0, 0, 0], [0, 0, 1]])) 
    ]
    for i, direction in enumerate(directions):
        fn.FunctionAssigner(sourceSpace.sub(i), sourceScalarSpace).assign(destinationSpace.sub(i), fn.project(fn.inner(sourceData, direction), sourceScalarSpace, solver_type='cg', preconditioner_type='hypre_amg'))
    
    return destinationSpace

Here are the results with the mumps solver (normal fenics.project()).

"""
Mesh size                                 Time to complete:
15x15x15                                               1.4s
25x25x25                                              33.3s
40x40x40                                             554.0s
55x55x55                                             Killed
"""

error output (myProjectFn())

Shapes do not match: <ComponentTensor id=140696883287600> and <Coefficient id=140696881005568>.
Traceback (most recent call last):
  File "testProjections.py", line 86, in <module>
    E = myProjectFn(sourceSpace, sourceScalarSpace, sourceData, destinationSpace)
  File "testProjections.py", line 77, in myProjectFn
    fn.FunctionAssigner(sourceSpace.sub(i), sourceScalarSpace).assign(destinationSpace.sub(i), fn.project(fn.inner(sourceData, direction), sourceScalarSpace, solver_type='cg', preconditioner_type='hypre_amg'))
  File "/usr/lib/python3/dist-packages/ufl/operators.py", line 158, in inner
    return Inner(a, b)
  File "/usr/lib/python3/dist-packages/ufl/tensoralgebra.py", line 147, in __new__
    error("Shapes do not match: %s and %s." % (ufl_err_str(a), ufl_err_str(b)))
  File "/usr/lib/python3/dist-packages/ufl/log.py", line 158, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Shapes do not match: <ComponentTensor id=140696883287600> and <Coefficient id=140696881005568>.

Full MWE:

# Ryan Budde 2021
import fenics as fn
import time
import numpy

tic = time.perf_counter() # start timer

# define left and right walls
def left(x, on_boundary):
    return fn.near(x[0], 0) 
def right(x, on_boundary):
    return fn.near(x[0], 1)

# define mesh
mesh = fn.UnitCubeMesh(15,15,15)
saveFileName = 'coarse'
# comment or uncomment the below to run finer mesh
# which greatly increases compute time
# and kills the program
mesh = fn.UnitCubeMesh(55,55,55)
saveFileName = 'fine'

# define space
V = fn.FunctionSpace(mesh, 'CG', 1)
dx = fn.Measure('dx', domain=mesh)

bc1 = fn.DirichletBC(V, fn.Constant(0),  left)
bc2 = fn.DirichletBC(V, fn.Constant(1), right)
bcs = [bc1, bc2]

# equation setup
u = fn.TrialFunction(V)
v = fn.TestFunction(V)
a = fn.dot(fn.grad(u), fn.grad(v)) * fn.dx
L = fn.Constant('0') * v * fn.dx

# explicit solver (apparently more efficient)
A = fn.assemble(a)
b = fn.assemble(L)
for bc in bcs:
   bc.apply(A,b)
solver = fn.KrylovSolver('gmres')
prm = solver.parameters
prm['absolute_tolerance'] = 1E-7
prm['relative_tolerance'] = 1E-4
prm['maximum_iterations'] = 10000
u = fn.Function(V)
U = u.vector()

print(f"Begin solve\t{time.perf_counter()-tic:0.4f}")
solver.solve(A, U, b)
print(f"Solve complete\t{time.perf_counter()-tic:0.4f}")

# Calculate field from projection of potential
E = -fn.grad(u)
print(f"Start projection\t{time.perf_counter()-tic:0.4f}")

# E = fn.project(E, solver_type="mumps")

# attempt at solving in each dimension
####################################################################
def myProjectFn(sourceSpace, sourceScalarSpace, sourceData, destinationSpace):
    directions = [
        fn.Constant(numpy.array([[1, 0, 0], [0, 0, 0], [0, 0, 0]])) ,
        fn.Constant(numpy.array([[0, 0, 0], [0, 1, 0], [0, 0, 0]])) ,
        fn.Constant(numpy.array([[0, 0, 0], [0, 0, 0], [0, 0, 1]])) 
    ]
    for i, direction in enumerate(directions):
        fn.FunctionAssigner(sourceSpace.sub(i), sourceScalarSpace).assign(destinationSpace.sub(i), fn.project(fn.inner(sourceData, direction), sourceScalarSpace, solver_type='cg', preconditioner_type='hypre_amg'))
    
    return destinationSpace
sourceSpace = fn.VectorFunctionSpace(mesh, 'CG', 1)
sourceScalarSpace = fn.FunctionSpace(mesh, 'CG', 1)
destinationSpace = fn.Function(sourceSpace)
sourceData = E
E = myProjectFn(sourceSpace, sourceScalarSpace, sourceData, destinationSpace)
####################################################################

print(f"Projection complete\t{time.perf_counter()-tic:0.4f}")

# export files if desired
# pFile = fn.File(saveFileName + 'testProject_potential.pvd')
# pFile << u
# eFile = fn.File(saveFileName + 'testProject_e_field.pvd')
# eFile << E

print(f"Data exported, end of file\t{time.perf_counter()-tic:0.4f}")

Thanks to all for assistance.

2

I found a solution for my specific case based on the posts here. This solution works to calculate a DG 0 gradient from a CG 1 source. I think there still may be ways to explicitly solve with the fenics.project() function is ways that are both efficient and general, but this workaround will be sufficient for now. With this method the 55x55x55 mesh which previously killed the program (at least 12 minutes) now runs in about 3 seconds.

The gradient is calculated as

Q = fn.VectorFunctionSpace(mesh, 'DG', 0)
p = fn.TrialFunction(Q)
q = fn.TestFunction(Q)
L = fn.inner(-fn.grad(u), q)*dx # u is previously calculated potential solution
M = fn.assemble(fn.inner(p,q)*dx)
b = fn.assemble(L)
grad_u0 = fn.Function(Q) # grad_u0 is the final answer
x0 = grad_u0.vector()
fn.solve(M, x0, b, 'cg', 'hypre_amg')

Complete code:

# Ryan Budde 2021
import fenics as fn
import time
import numpy

tic = time.perf_counter() # start timer

# define left and right walls
def left(x, on_boundary):
    return fn.near(x[0], 0) 
def right(x, on_boundary):
    return fn.near(x[0], 1)

# define mesh
mesh = fn.UnitCubeMesh(15,15,15)
saveFileName = 'coarse'
# uncomment the below to run finer mesh
# which greatly increases compute time
# and kills the program
mesh = fn.UnitCubeMesh(55,55,55)
saveFileName = 'fine'

# define space
V = fn.FunctionSpace(mesh, 'CG', 1)
dx = fn.Measure('dx', domain=mesh)

bc1 = fn.DirichletBC(V, fn.Constant(0),  left)
bc2 = fn.DirichletBC(V, fn.Constant(1), right)
bcs = [bc1, bc2]

# equation setup
u = fn.TrialFunction(V)
v = fn.TestFunction(V)
a = fn.dot(fn.grad(u), fn.grad(v)) * fn.dx
L = fn.Constant('0') * v * fn.dx

# explicit solver (apparently more efficient)
A = fn.assemble(a)
b = fn.assemble(L)
for bc in bcs:
   bc.apply(A,b)
solver = fn.KrylovSolver('gmres')
prm = solver.parameters
prm['absolute_tolerance'] = 1E-7
prm['relative_tolerance'] = 1E-4
prm['maximum_iterations'] = 10000
u = fn.Function(V)
U = u.vector()
solver.solve(A, U, b)


# solve gradient
print(f"Start gradient projection\t{time.perf_counter()-tic:0.4f}")
Q = fn.VectorFunctionSpace(mesh, 'DG', 0)
p = fn.TrialFunction(Q)
q = fn.TestFunction(Q)
L = fn.inner(-fn.grad(u), q)*dx
M = fn.assemble(fn.inner(p,q)*dx)
b = fn.assemble(L)
grad_u0 = fn.Function(Q)
x0 = grad_u0.vector()
fn.solve(M, x0, b, 'cg', 'hypre_amg')
print(f"Grad projection complete\t{time.perf_counter()-tic:0.4f}")

# alternative, which takes forever
# grad_u0 = fn.project(-fn.grad(u), solver_type="mumps")

# export files if desired
pFile = fn.File(saveFileName + 'testProject_potential.pvd')
pFile << u
eFile = fn.File(saveFileName + 'testProject_e_field.pvd')
eFile << grad_u0
2 Likes