Skip to content

Commit b15d632

Browse files
RemDelaporteMathurinRemDelaporteMathurin
committed
unified discontinuous classes
1 parent c072a5c commit b15d632

14 files changed

+96
-168
lines changed

‎examples/advection_diffusion_multi_material.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -97,7 +97,7 @@ def velocity_func(x):
9797
writer = dolfinx.io.VTXWriter(MPI.COMM_WORLD, "velocity.bp", u, engine="BP5")
9898
writer.write(t=0)
9999

100-
my_model = F.HTransportProblemDiscontinuous()
100+
my_model = F.HydrogenTransportProblemDiscontinuous()
101101
my_model.mesh = F.Mesh(mesh)
102102
my_model.volume_meshtags = ct
103103
my_model.facet_meshtags = mt

‎examples/multi_material_1d.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@
22

33
import festim as F
44

5-
my_model = F.HTransportProblemDiscontinuous()
5+
my_model = F.HydrogenTransportProblemDiscontinuous()
66

77
interface_1 = 0.5
88
interface_2 = 0.7

‎examples/multi_material_2d.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -63,7 +63,7 @@ def half(x):
6363

6464
mesh, mt, ct = generate_mesh()
6565

66-
my_model = F.HTransportProblemDiscontinuous()
66+
my_model = F.HydrogenTransportProblemDiscontinuous()
6767
my_model.mesh = F.Mesh(mesh)
6868
my_model.volume_meshtags = ct
6969
my_model.facet_meshtags = mt

‎examples/multi_material_2d_bis.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -63,7 +63,7 @@ def half(x):
6363

6464
mesh, mt, ct = generate_mesh()
6565

66-
my_model = F.HTransportProblemDiscontinuous()
66+
my_model = F.HydrogenTransportProblemDiscontinuous()
6767
my_model.mesh = F.Mesh(mesh)
6868
my_model.volume_meshtags = ct
6969
my_model.facet_meshtags = mt

‎examples/multi_material_transient.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@
22

33
import festim as F
44

5-
my_model = F.HTransportProblemDiscontinuous()
5+
my_model = F.HydrogenTransportProblemDiscontinuous()
66

77
interface_1 = 0.5
88
interface_2 = 0.7

‎examples/multi_material_with_one_material.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@
22

33
import festim as F
44

5-
my_model = F.HTransportProblemDiscontinuous()
5+
my_model = F.HydrogenTransportProblemDiscontinuous()
66

77

88
N = 1500

‎src/festim/__init__.py‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -47,8 +47,7 @@
4747
)
4848
from .problem import ProblemBase
4949
from .hydrogen_transport_problem import (
50-
HTransportProblemDiscontinuous,
51-
HTransportProblemPenalty,
50+
HydrogenTransportProblemDiscontinuous,
5251
HydrogenTransportProblem,
5352
HydrogenTransportProblemDiscontinuousChangeVar,
5453
)

‎src/festim/coupled_heat_hydrogen_problem.py‎

Lines changed: 4 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -4,8 +4,7 @@
44
from festim.heat_transfer_problem import HeatTransferProblem
55
from festim.helpers import as_fenics_constant, nmm_interpolate
66
from festim.hydrogen_transport_problem import (
7-
HTransportProblemDiscontinuous,
8-
HTransportProblemPenalty,
7+
HydrogenTransportProblemDiscontinuous,
98
HydrogenTransportProblem,
109
HydrogenTransportProblemDiscontinuousChangeVar,
1110
)
@@ -91,15 +90,13 @@ def hydrogen_problem(self):
9190
def hydrogen_problem(self, value):
9291
if isinstance(
9392
value,
94-
HTransportProblemDiscontinuous
95-
| HTransportProblemPenalty
93+
HydrogenTransportProblemDiscontinuous
9694
| HydrogenTransportProblemDiscontinuousChangeVar,
9795
):
9896
raise NotImplementedError(
9997
"Coupled heat transfer - hydrogen transport simulations with "
100-
"HydrogenTransportProblemDiscontinuousChangeVar, "
101-
"HTransportProblemPenalty or"
102-
"HydrogenTransportProblemDiscontinuousChangeVar, "
98+
"HydrogenTransportProblemDiscontinuousChangeVar or"
99+
"HydrogenTransportProblemDiscontinuous"
103100
"not currently supported"
104101
)
105102
elif not isinstance(value, HydrogenTransportProblem):

‎src/festim/hydrogen_transport_problem.py‎

Lines changed: 75 additions & 141 deletions
Original file line numberDiff line numberDiff line change
@@ -36,7 +36,7 @@
3636
from festim.helpers import as_fenics_constant, get_interpolation_points
3737
from festim.mesh import Mesh
3838

39-
__all__ = ["HTransportProblemDiscontinuous", "HydrogenTransportProblem"]
39+
__all__ = ["HydrogenTransportProblemDiscontinuous", "HydrogenTransportProblem"]
4040

4141

4242
class HydrogenTransportProblem(problem.ProblemBase):
@@ -948,10 +948,10 @@ def post_processing(self):
948948
vtxfile.write(float(self.t))
949949

950950

951-
class HTransportProblemDiscontinuous(HydrogenTransportProblem):
951+
class HydrogenTransportProblemDiscontinuous(HydrogenTransportProblem):
952952
interfaces: list[_subdomain.Interface]
953-
petsc_options: dict
954953
surface_to_volume: dict
954+
method_interface: str = "penalty"
955955

956956
def __init__(
957957
self,
@@ -1002,15 +1002,10 @@ def __init__(
10021002
settings,
10031003
exports,
10041004
traps,
1005+
petsc_options=petsc_options,
10051006
)
10061007
self.interfaces = interfaces or []
10071008
self.surface_to_volume = surface_to_volume or {}
1008-
default_petsc_options = {
1009-
"ksp_type": "preonly",
1010-
"pc_type": "lu",
1011-
"pc_factor_mat_solver_type": "mumps",
1012-
}
1013-
self.petsc_options = petsc_options or default_petsc_options
10141009
self._vtxfiles: list[dolfinx.io.VTXWriter] = []
10151010

10161011
def initialise(self):
@@ -1343,36 +1338,79 @@ def mixed_term(u, v, n):
13431338
self.mesh.mesh, self.temperature_fenics(res[1]), H
13441339
)
13451340

1346-
F_0 = -0.5 * mixed_term((u_b + u_t), v_b, n_0) * dInterface(
1347-
interface.id
1348-
) - 0.5 * mixed_term(v_b, (u_b / K_b - u_t / K_t), n_0) * dInterface(
1349-
interface.id
1350-
)
1341+
if self.method_interface == "penalty":
1342+
if (
1343+
subdomain_0.material.solubility_law
1344+
== subdomain_1.material.solubility_law
1345+
):
1346+
left = u_b / K_b
1347+
right = u_t / K_t
1348+
else:
1349+
if subdomain_0.material.solubility_law == "henry":
1350+
left = u_b / K_b
1351+
elif subdomain_0.material.solubility_law == "sievert":
1352+
left = (u_b / K_b) ** 2
1353+
else:
1354+
raise ValueError(
1355+
f"Unknown material law {subdomain_0.material.solubility_law}"
1356+
)
13511357

1352-
F_1 = +0.5 * mixed_term((u_b + u_t), v_t, n_0) * dInterface(
1353-
interface.id
1354-
) - 0.5 * mixed_term(v_t, (u_b / K_b - u_t / K_t), n_0) * dInterface(
1355-
interface.id
1356-
)
1357-
F_0 += (
1358-
2
1359-
* gamma
1360-
/ (h_0 + h_1)
1361-
* (u_b / K_b - u_t / K_t)
1362-
* v_b
1363-
* dInterface(interface.id)
1364-
)
1365-
F_1 += (
1366-
-2
1367-
* gamma
1368-
/ (h_0 + h_1)
1369-
* (u_b / K_b - u_t / K_t)
1370-
* v_t
1371-
* dInterface(interface.id)
1372-
)
1358+
if subdomain_1.material.solubility_law == "henry":
1359+
right = u_t / K_t
1360+
elif subdomain_1.material.solubility_law == "sievert":
1361+
right = (u_t / K_t) ** 2
1362+
else:
1363+
raise ValueError(
1364+
f"Unknown material law {subdomain_1.material.solubility_law}"
1365+
)
1366+
1367+
equality = right - left
13731368

1374-
subdomain_0.F += F_0
1375-
subdomain_1.F += F_1
1369+
F_0 = (
1370+
interface.penalty_term
1371+
* ufl.inner(equality, v_b)
1372+
* dInterface(interface.id)
1373+
)
1374+
F_1 = (
1375+
-interface.penalty_term
1376+
* ufl.inner(equality, v_t)
1377+
* dInterface(interface.id)
1378+
)
1379+
1380+
subdomain_0.F += F_0
1381+
subdomain_1.F += F_1
1382+
1383+
elif self.method_interface == "nietsche":
1384+
F_0 = -0.5 * mixed_term((u_b + u_t), v_b, n_0) * dInterface(
1385+
interface.id
1386+
) - 0.5 * mixed_term(v_b, (u_b / K_b - u_t / K_t), n_0) * dInterface(
1387+
interface.id
1388+
)
1389+
1390+
F_1 = +0.5 * mixed_term((u_b + u_t), v_t, n_0) * dInterface(
1391+
interface.id
1392+
) - 0.5 * mixed_term(v_t, (u_b / K_b - u_t / K_t), n_0) * dInterface(
1393+
interface.id
1394+
)
1395+
F_0 += (
1396+
2
1397+
* gamma
1398+
/ (h_0 + h_1)
1399+
* (u_b / K_b - u_t / K_t)
1400+
* v_b
1401+
* dInterface(interface.id)
1402+
)
1403+
F_1 += (
1404+
-2
1405+
* gamma
1406+
/ (h_0 + h_1)
1407+
* (u_b / K_b - u_t / K_t)
1408+
* v_t
1409+
* dInterface(interface.id)
1410+
)
1411+
1412+
subdomain_0.F += F_0
1413+
subdomain_1.F += F_1
13761414

13771415
J = []
13781416
# this is the symbolic differentiation of the Jacobian
@@ -1520,110 +1558,6 @@ def __del__(self):
15201558
vtxfile.close()
15211559

15221560

1523-
class HTransportProblemPenalty(HTransportProblemDiscontinuous):
1524-
def create_formulation(self):
1525-
"""
1526-
Takes all the formulations for each subdomain and adds the interface conditions.
1527-
1528-
Finally compute the jacobian matrix and store it in the ``J`` attribute,
1529-
adds the ``entity_maps`` to the forms and store them in the ``forms`` attribute
1530-
"""
1531-
mesh = self.mesh.mesh
1532-
mt = self.facet_meshtags
1533-
1534-
for interface in self.interfaces:
1535-
interface.mesh = mesh
1536-
interface.mt = mt
1537-
1538-
integral_data = [
1539-
interface.compute_mapped_interior_facet_data(mesh)
1540-
for interface in self.interfaces
1541-
]
1542-
[interface.pad_parent_maps() for interface in self.interfaces]
1543-
dInterface = ufl.Measure("dS", domain=mesh, subdomain_data=integral_data)
1544-
1545-
entity_maps = {
1546-
sd.submesh: sd.parent_to_submesh for sd in self.volume_subdomains
1547-
}
1548-
for interface in self.interfaces:
1549-
subdomain_0, subdomain_1 = interface.subdomains
1550-
res = interface.restriction
1551-
1552-
all_mobile_species = [spe for spe in self.species if spe.mobile]
1553-
if len(all_mobile_species) > 1:
1554-
raise NotImplementedError("Multiple mobile species not implemented")
1555-
H = all_mobile_species[0]
1556-
v_b = H.subdomain_to_test_function[subdomain_0](res[0])
1557-
v_t = H.subdomain_to_test_function[subdomain_1](res[1])
1558-
1559-
u_b = H.subdomain_to_solution[subdomain_0](res[0])
1560-
u_t = H.subdomain_to_solution[subdomain_1](res[1])
1561-
1562-
K_b = subdomain_0.material.get_solubility_coefficient(
1563-
self.mesh.mesh, self.temperature_fenics(res[0]), H
1564-
)
1565-
K_t = subdomain_1.material.get_solubility_coefficient(
1566-
self.mesh.mesh, self.temperature_fenics(res[1]), H
1567-
)
1568-
1569-
if (
1570-
subdomain_0.material.solubility_law
1571-
== subdomain_1.material.solubility_law
1572-
):
1573-
left = u_b / K_b
1574-
right = u_t / K_t
1575-
else:
1576-
if subdomain_0.material.solubility_law == "henry":
1577-
left = u_b / K_b
1578-
elif subdomain_0.material.solubility_law == "sievert":
1579-
left = (u_b / K_b) ** 2
1580-
else:
1581-
raise ValueError(
1582-
f"Unknown material law {subdomain_0.material.solubility_law}"
1583-
)
1584-
1585-
if subdomain_1.material.solubility_law == "henry":
1586-
right = u_t / K_t
1587-
elif subdomain_1.material.solubility_law == "sievert":
1588-
right = (u_t / K_t) ** 2
1589-
else:
1590-
raise ValueError(
1591-
f"Unknown material law {subdomain_1.material.solubility_law}"
1592-
)
1593-
1594-
equality = right - left
1595-
1596-
F_0 = (
1597-
interface.penalty_term
1598-
* ufl.inner(equality, v_b)
1599-
* dInterface(interface.id)
1600-
)
1601-
F_1 = (
1602-
-interface.penalty_term
1603-
* ufl.inner(equality, v_t)
1604-
* dInterface(interface.id)
1605-
)
1606-
1607-
subdomain_0.F += F_0
1608-
subdomain_1.F += F_1
1609-
1610-
J = []
1611-
# this is the symbolic differentiation of the Jacobian
1612-
for subdomain1 in self.volume_subdomains:
1613-
jac = []
1614-
for subdomain2 in self.volume_subdomains:
1615-
jac.append(
1616-
ufl.derivative(subdomain1.F, subdomain2.u),
1617-
)
1618-
J.append(jac)
1619-
# compile jacobian (J) and residual (F)
1620-
self.forms = dolfinx.fem.form(
1621-
[subdomain.F for subdomain in self.volume_subdomains],
1622-
entity_maps=entity_maps,
1623-
)
1624-
self.J = dolfinx.fem.form(J, entity_maps=entity_maps)
1625-
1626-
16271561
class HydrogenTransportProblemDiscontinuousChangeVar(HydrogenTransportProblem):
16281562
species: List[_species.Species]
16291563

‎src/festim/subdomain/volume_subdomain.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -60,7 +60,7 @@ def create_subdomain(self, mesh: dolfinx.mesh.Mesh, marker: dolfinx.mesh.MeshTag
6060
Creates the following attributes: ``.parent_mesh``, ``.submesh``, ``.submesh_to_mesh``,
6161
``.v_map``, ``padded``, and the entity map ``parent_to_submesh``.
6262
63-
Only used in ``festim.HTransportProblemDiscontinuous``
63+
Only used in ``festim.HydrogenTransportProblemDiscontinuous``
6464
6565
Args:
6666
mesh (dolfinx.mesh.Mesh): the parent mesh

0 commit comments

Comments
 (0)