[RemDelaporteMathurin]

A homogenised polycrystal, a rolled or columnar microstructure, and any non-cubic lattice conduct differently along different directions, and the property that describes that is a second-rank tensor. D_0 may now be given as a square matrix:

F.Material(D_0=[[9e-11, 2e-12], [2e-12, 2.8e-11]], E_D=0.1)

E_D stays a scalar -- one activation energy shared by every direction, the prefactor carrying the anisotropy. For per-direction activation energies, which a matrix D_0 deliberately cannot express, Material.D now also accepts a fem.Constant, a ufl expression or a plain matrix as well as a fem.Function, so a tensor built from the problem's temperature stays temperature dependent.

Four places read D as a scalar and had to be written with the matrix product. Each reduces to the previous expression when D is a scalar:

• SurfaceFlux: -dot(D grad(u), n) rather than -D dot(grad(u), n), which for a tensor cannot even be assembled.
• Nitsche: the consistency and symmetry terms take the matrix product, and the penalty scales on the normal conductance n.D.n, which is the conductance the boundary actually sees.
• define_D_global builds tensor-valued DG spaces when the material is anisotropic.
• as_fenics_constant builds tensor constants.

ChangeVar and the drift terms already used D * grad(...) and are unchanged.

Also fixes a latent bug this feature would otherwise trip over: if self.D_0 and self.D in Material.__init__ raises "truth value of an array is ambiguous" for any array input, as does if material.D in define_D_global. Both now test against None.

Tests: a manufactured solution sin(pi x) sin(pi y) on the unit square, run with a rotated tensor as well as a diagonal one -- with D_xy == 0 the off-diagonal term of the source vanishes and a formulation that mishandled it would still pass. Second-order convergence for identity, diagonal, rotated and strongly rotated tensors, the Arrhenius factor applied to the tensor, the surface flux of a slab against its analytical value and against the same integral assembled by hand, and weak Dirichlet enforcement against strong.

Description

Summary

Related Issues

Closes #754

Motivation and Context

Type of Change

• 🐛 Bug fix (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 (no functional changes, no API changes)
• 📝 Documentation update
• ✅ Test update (adding missing tests or correcting existing tests)
• 🔧 Build/CI configuration change

Testing

• All existing tests pass locally (pytest)
• I have added new tests that prove my fix is effective or that my feature works

Code Quality Checklist

• My code follows the code style of this project (Ruff formatted: ruff format .)
• My code passes linting checks (ruff check .)
• I have performed a self-review of my own code
• I have commented my code, particularly in hard-to-understand areas

Documentation

• I have updated the documentation accordingly (if applicable)
• I have added docstrings to new functions/classes following the project conventions

Breaking Changes

Screenshots/Examples

Additional Notes

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 96.11%. Comparing base (ad055a8) to head (7112000).

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1260      +/-   ##
==========================================
+ Coverage   96.07%   96.11%   +0.04%     
==========================================
  Files          57       57              
  Lines        5069     5126      +57     
==========================================
+ Hits         4870     4927      +57     
  Misses        199      199              

Flag	Coverage Δ
dolfinx-v0.10.0	93.21% <100.00%> (+0.07%)	⬆️
dolfinx-v0.11.0	95.98% <100.00%> (+0.04%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Harness.
📢 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.

[ee-nn]

@RemDelaporteMathurin Let me know if you want to finish this at some point, otherwise I could take over!

[RemDelaporteMathurin]

Sure if you want to you can take this over! Thanks!

[ee-nn]

I feel like the sooner #1191 can be implemented the better because all of this is going to be axed sooner or later anyway. Perhaps I can provide a PR for that soon

[ee-nn]

@RemDelaporteMathurin I think this should more or less be good now, let me know if you have any comments or changes! I tested out the new interface with a small mwe and also I can confirm that all the tests pass

from mpi4py import MPI

import dolfinx
import numpy as np

import festim as F

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 32, 32)

# D_0 is a tensor in the mesh coordinate system: R diag(D_fast, D_slow) R^T.
angle = np.deg2rad(30.0)
rotation = np.array([[np.cos(angle), -np.sin(angle)], [np.sin(angle), np.cos(angle)]])
D_0 = rotation @ np.diag([5e-7, 5e-8]) @ rotation.T
E_D = 0.2
T = 500.0
material = F.Material(D_0=D_0, E_D=E_D)
# Alternatively, one could do F.Material(D=D_0 * np.exp(-E_D / (F.k_B * T)))

volume = F.VolumeSubdomain(
    id=1,
    material=material,
    locator=lambda x: np.ones_like(x[0], dtype=bool),
)
left = F.SurfaceSubdomain(id=2, locator=lambda x: np.isclose(x[0], 0.0))
right = F.SurfaceSubdomain(id=3, locator=lambda x: np.isclose(x[0], 1.0))
hydrogen = F.Species("H", subdomains=[volume])
inlet_flux = F.SurfaceFlux(field=hydrogen, surface=left)
outlet_flux = F.SurfaceFlux(field=hydrogen, surface=right)

model = F.HydrogenTransportProblemDiscontinuous(
    mesh=F.Mesh(mesh),
    species=[hydrogen],
    subdomains=[volume, left, right],
    boundary_conditions=[
        F.FixedConcentrationBC(subdomain=left, value=1.0, species=hydrogen),
        F.FixedConcentrationBC(subdomain=right, value=0.0, species=hydrogen),
    ],
    temperature=T,
    settings=F.Settings(transient=False, atol=1e-14, rtol=1e-10),
    exports=[inlet_flux, outlet_flux],
)
model.initialise()
model.run()

print(f"D at {T:g} K (m^2/s):\n{D_0 * np.exp(-E_D / (F.k_B * T))}")
# SurfaceFlux integrates -(D grad(c)) . n; positive means leaving.
print(f"Integrated inlet flux:  {inlet_flux.value:.6e}")
print(f"Integrated outlet flux: {outlet_flux.value:.6e}")