Surface Flux fix

[jhdark]

Proposed changes

With this change, the mesh normal needed to evaluate the surface flux can be evaluated from the field rather than passing it in Hydrogen or HeatTransferProblem. The method is different depending if the problem being solved is multispecies or not, as the solution will either be a function or an index.

This should make it easier to process each future surface quantity export as each one will only need the ds ufl.measure:

for export in self.exports:
    if isinstance(export, F.SurfaceQuantity):
        export.compute(self.ds)

Types of changes

What types of changes does your code introduce to FESTIM?

• Bugfix (non-breaking change which fixes an issue)
• New feature (non-breaking change which adds functionality)
• Breaking change (fix or feature that would cause existing functionality to not work as expected)
• Code refactoring
• Documentation Update (if none of the other choices apply)
• New tests

Checklist

• Black formatted
• Unit tests pass locally with my changes
• I have added tests that prove my fix is effective or that my feature works
• I have added necessary documentation (if appropriate)

Further comments

If this is a relatively large or complex change, kick off the discussion by explaining why you chose the solution you did and what alternatives you considered, etc...

[jhdark]

Not sure why the conda CI is failing, I can re-run the tests tomorrow and then it should be grand

[RemDelaporteMathurin]

This would work. However, FacetNormal would be computed at each write of each derived quantity.

For large meshes, computing it is expensive. We could test this with a very large 3D mesh, and do something like:

mesh = large_3d_mesh
u = Function()

# time this loop
for i in range(1000):
    n = ufl.FacetNormal(mesh)
    assemble_scalar(fem.form(ufl.dot(ufl.grad(self.field.solution), n) * ds))

as opposed to

mesh = large_3d_mesh
u = Function()
n = ufl.FacetNormal(mesh)
# time this loop
for i in range(1000):
    assemble_scalar(fem.form(ufl.dot(ufl.grad(self.field.solution), n) * ds))

I think that, since it is unique to the mesh, it ought to stay in mesh and is simply passed to compute

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 98.43%. Comparing base (aac14ac) to head (92f5798).

Additional details and impacted files

@@             Coverage Diff             @@
##           fenicsx     #747      +/-   ##
===========================================
- Coverage    98.49%   98.43%   -0.07%     
===========================================
  Files           28       28              
  Lines         1532     1536       +4     
===========================================
+ Hits          1509     1512       +3     
- Misses          23       24       +1     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[RemDelaporteMathurin]

@jhdark did some timing experiments and computing FacetNormal doesn't seem expensive at all even on big meshes

[jhdark]

For reference these were the test files I used:

from dolfinx.mesh import create_unit_cube, meshtags
from mpi4py import MPI
import ufl
from dolfinx import fem
import numpy as np
import time

test_mesh = create_unit_cube(MPI.COMM_WORLD, 100, 100, 200)
x = ufl.SpatialCoordinate(test_mesh)
facet_indices = np.array([], dtype=np.int32)
facet_tags = np.array([], dtype=np.int32)
facet_meshtags = meshtags(test_mesh, 2, facet_indices, facet_tags)
ds = ufl.Measure("ds", domain=test_mesh, subdomain_data=facet_meshtags)
V = fem.FunctionSpace(test_mesh, ("CG", 1))
u = fem.Function(V)

# time this loop
start_time = time.time()
for i in range(1000):
    n = ufl.FacetNormal(test_mesh)
    fem.assemble_scalar(fem.form(ufl.dot(ufl.grad(u), n) * ds))

time_with = time.time() - start_time

print("Time with: ", time_with)

and:

from dolfinx.mesh import create_unit_cube, meshtags
from mpi4py import MPI
import ufl
from dolfinx import fem
import numpy as np
import time

test_mesh = create_unit_cube(MPI.COMM_WORLD, 100, 100, 200)
x = ufl.SpatialCoordinate(test_mesh)
facet_indices = np.array([], dtype=np.int32)
facet_tags = np.array([], dtype=np.int32)
facet_meshtags = meshtags(test_mesh, 2, facet_indices, facet_tags)
ds = ufl.Measure("ds", domain=test_mesh, subdomain_data=facet_meshtags)
V = fem.FunctionSpace(test_mesh, ("CG", 1))
u = fem.Function(V)
n = ufl.FacetNormal(test_mesh)

# time this loop
start_time_1 = time.time()
for i in range(1000):
    fem.assemble_scalar(fem.form(ufl.dot(ufl.grad(u), n) * ds))

time_without = time.time() - start_time_1

print("Time without: ", time_without)