Advection term support

[jhdark]

Proposed changes

Describe the big picture of your changes here to communicate to the maintainers why we should accept this pull request. If it fixes a bug or resolves a feature request, be sure to link to that issue.

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

Example usage in MMS test:

test_mesh_2d = create_unit_square(MPI.COMM_WORLD, 200, 200)
    x_2d = ufl.SpatialCoordinate(test_mesh_2d)

    # coupled simulation properties
    D_0, E_D = 1.2, 0.1
    k_0, E_k = 2.2, 0.6
    p_0, E_p = 0.5, 0.1
    n_trap = 100
    k_B = F.k_B

    # common festim objects
    test_mat = F.Material(D_0=D_0, E_D=E_D)
    test_vol_sub = F.VolumeSubdomain(id=1, material=test_mat)

    boundary = F.SurfaceSubdomain(id=1)
    test_mobile = F.Species("mobile", mobile=True)
    test_trapped = F.Species(name="trapped", mobile=True)
    empty_trap = F.ImplicitSpecies(n=n_trap, others=[test_trapped])

    V = fem.functionspace(test_mesh_2d, ("Lagrange", 1))
    T = fem.Function(V)
    T_expr = lambda x: 100 + 200 * x[0] + 100 * x[1]
    T.interpolate(T_expr)

    # create velocity field
    v_cg = basix.ufl.element(
        "Lagrange",
        test_mesh_2d.topology.cell_name(),
        2,
        shape=(test_mesh_2d.geometry.dim,),
    )
    V_velocity = fem.functionspace(test_mesh_2d, v_cg)
    u = fem.Function(V_velocity)

    def velocity_func(x):
        values = np.zeros((2, x.shape[1]))  # Initialize with zeros

        scalar_value = x[1] * (x[1] - 1)  # Compute the scalar function
        values[0] = scalar_value  # Assign to first component
        values[1] = 0  # Second component remains zero

        return values

    u.interpolate(velocity_func)

    # define hydrogen problem
    exact_mobile_solution = lambda x: 200 * x[0] ** 2 + 300 * x[1] ** 2
    exact_trapped_solution = lambda x: 10 * x[0] ** 2 + 10 * x[1] ** 2

    D = D_0 * ufl.exp(-E_D / (k_B * T))
    k = k_0 * ufl.exp(-E_k / (k_B * T))
    p = p_0 * ufl.exp(-E_p / (k_B * T))

    f = (
        -ufl.div(D * ufl.grad(exact_mobile_solution(x_2d)))
        + ufl.inner(u, ufl.grad(exact_mobile_solution(x_2d)))
        + k * exact_mobile_solution(x_2d) * (n_trap - exact_trapped_solution(x_2d))
        - p * exact_trapped_solution(x_2d)
    )

    g = (
        -ufl.div(D * ufl.grad(exact_trapped_solution(x_2d)))
        - k * exact_mobile_solution(x_2d) * (n_trap - exact_trapped_solution(x_2d))
        + p * exact_trapped_solution(x_2d)
    )

    my_bcs = []
    for species, value in zip(
        [test_mobile, test_trapped], [exact_mobile_solution, exact_trapped_solution]
    ):
        my_bcs.append(
            F.FixedConcentrationBC(subdomain=boundary, value=value, species=species)
        )

    test_hydrogen_problem = F.HydrogenTransportProblem(
        mesh=F.Mesh(test_mesh_2d),
        subdomains=[test_vol_sub, boundary],
        boundary_conditions=my_bcs,
        species=[test_mobile, test_trapped],
        temperature=T,
        reactions=[
            F.Reaction(
                reactant=[test_mobile, empty_trap],
                product=test_trapped,
                k_0=k_0,
                E_k=E_k,
                p_0=p_0,
                E_p=E_p,
                volume=test_vol_sub,
            )
        ],
        sources=[
            F.ParticleSource(value=f, volume=test_vol_sub, species=test_mobile),
            F.ParticleSource(value=g, volume=test_vol_sub, species=test_trapped),
        ],
        advection_terms=[
            F.AdvectionTerm(velocity=u, subdomain=test_vol_sub, species=test_mobile)
        ],
        settings=F.Settings(
            atol=1e-10,
            rtol=1e-10,
            transient=False,
        ),
    )

    test_hydrogen_problem.initialise()
    test_hydrogen_problem.run()

    # compare computed values with exact solutions
    mobile_computed = test_mobile.post_processing_solution
    trapped_computed = test_trapped.post_processing_solution

    L2_error_mobile = error_L2(mobile_computed, exact_mobile_solution)
    L2_error_trapped = error_L2(trapped_computed, exact_trapped_solution)

[codecov]

Codecov Report

Attention: Patch coverage is 97.29730% with 3 lines in your changes missing coverage. Please review.

Project coverage is 96.14%. Comparing base (dcf7a01) to head (cf5232e).
Report is 191 commits behind head on fenicsx.

Files with missing lines	Patch %	Lines
src/festim/hydrogen_transport_problem.py	91.42%	3 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##           fenicsx     #965      +/-   ##
===========================================
+ Coverage    96.09%   96.14%   +0.04%     
===========================================
  Files           46       47       +1     
  Lines         2611     2721     +110     
===========================================
+ Hits          2509     2616     +107     
- Misses         102      105       +3     

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

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[RemDelaporteMathurin]

Looking good! Just wondering if we could test more the multi-material case

[RemDelaporteMathurin]

I know we want to merge this but should we take the time to actually test something here?

[jhdark]

Yeah we could have a proper MMS test for this one too

[RemDelaporteMathurin]

We should at least write it down as an issue to remind us of it

[jhdark]

Made a note that we should add another MMS test in a discontinuous case in #969