Efficient way to compute surface gradient

Question

Dear Fellow-Fenics-Users,

I am using Fenics/Dolfin to solve a steady state advection diffusion equation.
I am eventually interested in the concentration gradient at (a part of the) domain boundary.

At the moment, I compute this gradient as follows

(1) Compute gradient of c in whole domain:

V_g = VectorFunctionSpace(mesh, 'Lagrange', 1)
w = TrialFunction(V_g)
v = TestFunction(V_g)

a = inner(w, v)*dx
L = inner(-grad(c), v)*dx
gradc = Function(V_g)
solve(a == L, gradc,solver_parameters={'linear_solver':'mumps'})

Note: I decided to use this form instead of

gradc = project(-grad(c),V_g)

to avoid problems with memory usage in large domains.

(2) Compute the integral over the relevant surface:

dA = ds(subdomain_data=Boundaries)
Area = assemble(1*dA(2,domain=mesh))
TotalFlux = assemble(gradc[0]*dA(2,domain=mesh))
Average Flux = TotalFlux/Area

Note that the surface in question lays in the (y,z) plane, so the x-component of the gradient is the flux I am interested in.

However, computing the gradient this way is incredibly slow. Is there anyway to avoid computing the gradient for the whole domain and to consider only the gradient at the surface?

Thank you very much in advance and best wishes,
Merlin :)

Answer 1

Hi, I haven't seen a way of finding the area, so assembling 1 sounds fine. For the flux, if you already have c, have you tried integrating with Gauss' theorem? Because you could do

$$ \int_{\partial \Omega} \nabla c \cdot \vec e_1 = \int_\Omega \nabla\cdot (\nabla c\cdot \vec e_1) $$

and so try something like

assemble( div( c.dx(0))*dx, mesh = relevant_mesh ) 

Hope it works.