Skip to content

Fix penalty term #1030

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 8 commits into
base: fenicsx
Choose a base branch
from
Open

Fix penalty term #1030

Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
be376ac
started adding new penalty
RemDelaporteMathurin Sep 17, 2025
9537ccc
working now!
RemDelaporteMathurin Sep 17, 2025
aaff7f2
fix subdomain_a
RemDelaporteMathurin Sep 17, 2025
a2b6255
Merge branch 'fenicsx' into fix-penalty
RemDelaporteMathurin Sep 17, 2025
b374648
untraced to unrestrcited
RemDelaporteMathurin Sep 17, 2025
b3da603
fixed nitsche
RemDelaporteMathurin Sep 17, 2025
956518d
Merge branch 'fix-penalty' of https://github.com/festim-dev/FESTIM in…
RemDelaporteMathurin Sep 17, 2025
a2958b9
Merge branch 'fenicsx' into fix-penalty
RemDelaporteMathurin Sep 19, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
132 changes: 90 additions & 42 deletions src/festim/hydrogen_transport_problem.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1521,91 +1521,139 @@ def mixed_term(u, v, n):
# TODO only do this if the species in defined in both domains of the
# interface?
for H in all_mobile_species:
v_b = H.subdomain_to_test_function[subdomain_0](res[0])
v_t = H.subdomain_to_test_function[subdomain_1](res[1])
v_0 = H.subdomain_to_test_function[subdomain_0](res[0])
v_1 = H.subdomain_to_test_function[subdomain_1](res[1])

u_b = H.subdomain_to_solution[subdomain_0](res[0])
u_t = H.subdomain_to_solution[subdomain_1](res[1])
u_0 = H.subdomain_to_solution[subdomain_0](res[0])
u_1 = H.subdomain_to_solution[subdomain_1](res[1])

K_b = subdomain_0.material.get_solubility_coefficient(
u_0_unrestricted = H.subdomain_to_solution[subdomain_0].ufl_operands[0]
u_1_unrestricted = H.subdomain_to_solution[subdomain_1].ufl_operands[0]

K_0 = subdomain_0.material.get_solubility_coefficient(
self.mesh.mesh, self.temperature_fenics(res[0]), H
)
K_t = subdomain_1.material.get_solubility_coefficient(
K_1 = subdomain_1.material.get_solubility_coefficient(
self.mesh.mesh, self.temperature_fenics(res[1]), H
)

match self.method_interface:
case InterfaceMethod.penalty:
# Add penalty term for alpha/2 * equation**2 dGamma
# omega_A is Sieverts
# omeaga_B is Henry

if (
subdomain_0.material.solubility_law
== subdomain_1.material.solubility_law
!= subdomain_1.material.solubility_law
):
left = u_b / K_b
right = u_t / K_t
else:
# find sieverts then use it for A
# use other one for B
match subdomain_0.material.solubility_law:
case SolubilityLaw.HENRY:
left = u_b / K_b
case SolubilityLaw.SIEVERT:
left = (u_b / K_b) ** 2
case _:
raise ValueError(
"Unsupported material law "
+ f"{subdomain_0.material.solubility_law}"
u_a, u_a_unrestricted, v_a, K_a = (
u_0,
u_0_unrestricted,
v_0,
K_0,
)

match subdomain_1.material.solubility_law:
u_b, u_b_unrestricted, v_b, K_b = (
u_1,
u_1_unrestricted,
v_1,
K_1,
)
subdomain_a = subdomain_0
subdomain_b = subdomain_1
case SolubilityLaw.HENRY:
right = u_t / K_t
case SolubilityLaw.SIEVERT:
right = (u_t / K_t) ** 2
case _:
raise ValueError(
f"Unsupported material law "
f"{subdomain_1.material.solubility_law}"
u_a, u_a_unrestricted, v_a, K_a = (
u_1,
u_1_unrestricted,
v_1,
K_1,
)
u_b, u_b_unrestricted, v_b, K_b = (
u_0,
u_0_unrestricted,
v_0,
K_0,
)
subdomain_a = subdomain_1
subdomain_b = subdomain_0

equality = right - left
n_sorption = 0.5
else:
u_a, u_a_unrestricted, v_a, K_a = (
u_0,
u_0_unrestricted,
v_0,
K_0,
)
u_b, u_b_unrestricted, v_b, K_b = (
u_1,
u_1_unrestricted,
v_1,
K_1,
)
n_sorption = 1
subdomain_a = subdomain_0
subdomain_b = subdomain_1

tol = dolfinx.fem.Constant(
mesh, 1e2 * np.finfo(dolfinx.default_scalar_type).eps
)

F_0 = (
cond_b = ufl.gt(abs(u_b), tol)
u_b_padded = ufl.conditional(cond_b, u_b, tol)
equation = u_a - K_a * (u_b_padded / K_b) ** n_sorption

deq_du_a = ufl.inner(
ufl.diff(equation, u_a_unrestricted)[0], v_a
)
deq_du_b = ufl.inner(
ufl.diff(equation, u_b_unrestricted)[0], v_b
)
F_a = (
interface.penalty_term
* ufl.inner(equality, v_b)
* equation
* deq_du_a
* dInterface(interface.id)
)
F_1 = (
-interface.penalty_term
* ufl.inner(equality, v_t)
F_b = (
interface.penalty_term
* equation
* deq_du_b
* dInterface(interface.id)
)

subdomain_0.F += F_0
subdomain_1.F += F_1
subdomain_a.F += F_a
subdomain_b.F += F_b
case InterfaceMethod.nitsche:
F_0 = -0.5 * mixed_term((u_b + u_t), v_b, n_0) * dInterface(
F_0 = -0.5 * mixed_term((u_0 + u_1), v_0, n_0) * dInterface(
interface.id
) - 0.5 * mixed_term(
v_b, (u_b / K_b - u_t / K_t), n_0
v_0, (u_0 / K_0 - u_1 / K_1), n_0
) * dInterface(interface.id)

F_1 = +0.5 * mixed_term((u_b + u_t), v_t, n_0) * dInterface(
F_1 = +0.5 * mixed_term((u_0 + u_1), v_1, n_0) * dInterface(
interface.id
) - 0.5 * mixed_term(
v_t, (u_b / K_b - u_t / K_t), n_0
v_1, (u_0 / K_0 - u_1 / K_1), n_0
) * dInterface(interface.id)
F_0 += (
2
* gamma
/ (h_0 + h_1)
* (u_b / K_b - u_t / K_t)
* v_b
* (u_0 / K_0 - u_1 / K_1)
* v_0
* dInterface(interface.id)
)
F_1 += (
-2
* gamma
/ (h_0 + h_1)
* (u_b / K_b - u_t / K_t)
* v_t
* (u_0 / K_0 - u_1 / K_1)
* v_1
* dInterface(interface.id)
)

Expand Down
Loading