Anisotropic diffusion coefficient - #1260
Draft
RemDelaporteMathurin wants to merge 2 commits into
Draft
Conversation
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.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #1260 +/- ##
==========================================
- Coverage 95.89% 95.86% -0.03%
==========================================
Files 57 57
Lines 5017 5058 +41
==========================================
+ Hits 4811 4849 +38
- Misses 206 209 +3
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
ee-nn
reviewed
Sep 12, 2026
| ) -> fem.Constant: | ||
| """Converts a value to a dolfinx.Constant. | ||
|
|
||
| Array-like values become tensor-valued constants, which is how an anisotropic |
Collaborator
There was a problem hiding this comment.
Suggested change
| Array-like values become tensor-valued constants, which is how an anisotropic | |
| Array-like values become tensor-valued `fem.Constant` objects, which is how an anisotropic |
Comment on lines
+37
to
50
| if isinstance(value, bool): | ||
| raise TypeError(f"Value must not be a bool, not {type(value)}") | ||
| if isinstance(value, float | int): | ||
| return fem.Constant(mesh, dolfinx.default_scalar_type(float(value))) | ||
| elif isinstance(value, fem.Constant): | ||
| return value | ||
| elif isinstance(value, np.ndarray | list | tuple): | ||
| array = np.asarray(value, dtype=dolfinx.default_scalar_type) | ||
| return fem.Constant(mesh, array) | ||
| else: | ||
| raise TypeError( | ||
| f"Value must be a float, an int or a dolfinx.Constant, not {type(value)}" | ||
| "Value must be a float, an int, an array-like or a dolfinx.Constant, " | ||
| f"not {type(value)}" | ||
| ) |
Collaborator
There was a problem hiding this comment.
The first isinstance check appears redundant with the else statement at the end of the conditional block
The TypeError statement itself is also confusing (...not be a bool, not type()...). If kept, it should be changed
Comment on lines
+152
to
+153
| "D must be a fem.Function, a fem.Constant, a ufl expression or an " | ||
| f"array-like (for an anisotropic tensor), not {type(value)}" |
Collaborator
There was a problem hiding this comment.
Suggested change
| "D must be a fem.Function, a fem.Constant, a ufl expression or an " | |
| f"array-like (for an anisotropic tensor), not {type(value)}" | |
| "D must be a fem.Function, fem.Constant, ufl expression or a " | |
| f"square array-like (for an anisotropic tensor), not {type(value)}" |
Comment on lines
+175
to
+181
| def is_anisotropic(self, species=None) -> bool: | ||
| """Whether this material's diffusivity is a tensor rather than a scalar. | ||
|
|
||
| Anything reading ``D`` as a scalar -- a surface flux, a Nitsche penalty -- | ||
| has to ask, because for a tensor the same expression has to be written | ||
| with the matrix product instead. | ||
| """ |
Collaborator
There was a problem hiding this comment.
Just because a diffusion coefficient is passed as a tensor doesn't necessarily mean that it's anisotropic, e.g.
|
|
||
| A homogenised polycrystal, a rolled or columnar microstructure, or any hcp | ||
| lattice conducts differently along different directions, and the material | ||
| property that describes that is a second-rank tensor. ``D_0`` may therefore be |
Collaborator
There was a problem hiding this comment.
In practice allowing the energy to vary with direction could be helpful, for example the energy to migrate across a defect could be different from the ideal lattice.
One could get away with factoring this into D_0 but only at constant temperature, I believe the workaround to this is to allow the user to pass an array of ufl expressions
| F.Material() | ||
|
|
||
|
|
||
| class TestAnisotropicDiffusion: |
Collaborator
There was a problem hiding this comment.
Is there a reason these are in their own class?
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
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_0may now be given as a square matrix:E_Dstays a scalar -- one activation energy shared by every direction, the prefactor carrying the anisotropy. For per-direction activation energies, which a matrixD_0deliberately cannot express,Material.Dnow also accepts afem.Constant, a ufl expression or a plain matrix as well as afem.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.n.D.n, which is the conductance the boundary actually sees.define_D_globalbuilds tensor-valued DG spaces when the material is anisotropic.as_fenics_constantbuilds tensor constants.ChangeVarand the drift terms already usedD * grad(...)and are unchanged.Also fixes a latent bug this feature would otherwise trip over:
if self.D_0 and self.DinMaterial.__init__raises "truth value of an array is ambiguous" for any array input, as doesif material.Dindefine_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 -- withD_xy == 0the 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
Motivation and Context
Type of Change
Testing
pytest)Code Quality Checklist
ruff format .)ruff check .)Documentation
Breaking Changes
Screenshots/Examples
Additional Notes