From be376acc59998d8efc9fa3641b3da5ae58dc230d Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 17 Sep 2025 14:23:38 -0400 Subject: [PATCH 1/5] started adding new penalty --- src/festim/hydrogen_transport_problem.py | 73 ++++++++++++++---------- 1 file changed, 42 insertions(+), 31 deletions(-) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index ed69110a9..16fbe6f2b 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -1521,60 +1521,71 @@ 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( + 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: + + def tr_a(v): + return v("+") + + def tr_b(v): + return v("-") + + # 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}" - ) - - match subdomain_1.material.solubility_law: + u_a, v_a, K_a = u_0, v_0, K_0 + u_b, v_b, K_b = u_1, v_1, K_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, v_a, K_a = u_1, v_1, K_1 + u_b, v_b, K_b = u_0, v_0, K_0 - equality = right - left + n_sorption = 0.5 + else: + u_a, v_a, K_a = u_0, v_0, K_0 + u_b, v_b, K_b = u_1, v_1, K_1 + n_sorption = 1 + + tol = dolfinx.fem.Constant( + mesh, 1e2 * np.finfo(dolfinx.default_scalar_type).eps + ) + cond_b = ufl.gt(abs(tr_b(u_b)), tol) + u_b_padded = ufl.conditional(cond_b, tr_b(u_b), tol) + equation = tr_a(u_a) - K_a * (u_b_padded / K_b) ** n_sorption + deq_du_a = ufl.diff(equation, u_a) * tr_a(v_a) + deq_du_b = ufl.diff(equation, u_b) * tr_b(v_b) F_0 = ( 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) + * equation + * deq_du_b * dInterface(interface.id) ) From 9537ccc519b624682f746e34a7fa47aaae032e49 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 17 Sep 2025 14:55:54 -0400 Subject: [PATCH 2/5] working now! --- src/festim/hydrogen_transport_problem.py | 93 ++++++++++++++++-------- 1 file changed, 62 insertions(+), 31 deletions(-) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 16fbe6f2b..2dab95178 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -1527,6 +1527,9 @@ def mixed_term(u, v, n): u_0 = H.subdomain_to_solution[subdomain_0](res[0]) u_1 = H.subdomain_to_solution[subdomain_1](res[1]) + 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 ) @@ -1536,13 +1539,6 @@ def mixed_term(u, v, n): match self.method_interface: case InterfaceMethod.penalty: - - def tr_a(v): - return v("+") - - def tr_b(v): - return v("-") - # Add penalty term for alpha/2 * equation**2 dGamma # omega_A is Sieverts # omeaga_B is Henry @@ -1555,68 +1551,103 @@ def tr_b(v): # use other one for B match subdomain_0.material.solubility_law: case SolubilityLaw.SIEVERT: - u_a, v_a, K_a = u_0, v_0, K_0 - u_b, v_b, K_b = u_1, v_1, K_1 + u_a, u_a_untraced, v_a, K_a = ( + u_0, + u_0_unrestricted, + v_0, + K_0, + ) + u_b, u_b_untraced, v_b, K_b = ( + u_1, + u_1_unrestricted, + v_1, + K_1, + ) + subdomain_a = subdomain_0 + subdomain_b = subdomain_1 case SolubilityLaw.HENRY: - u_a, v_a, K_a = u_1, v_1, K_1 - u_b, v_b, K_b = u_0, v_0, K_0 + u_a, u_a_untraced, v_a, K_a = ( + u_1, + u_1_unrestricted, + v_1, + K_1, + ) + u_b, u_b_untraced, v_b, K_b = ( + u_0, + u_0_unrestricted, + v_0, + K_0, + ) + subdomain_a = subdomain_1 + subdomain_b = subdomain_0 n_sorption = 0.5 else: - u_a, v_a, K_a = u_0, v_0, K_0 - u_b, v_b, K_b = u_1, v_1, K_1 + u_a, u_a_untraced, v_a, K_a = ( + u_0, + u_0_unrestricted, + v_0, + K_0, + ) + u_b, u_b_untraced, v_b, K_b = ( + u_1, + u_1_unrestricted, + v_1, + K_1, + ) n_sorption = 1 tol = dolfinx.fem.Constant( mesh, 1e2 * np.finfo(dolfinx.default_scalar_type).eps ) - cond_b = ufl.gt(abs(tr_b(u_b)), tol) - u_b_padded = ufl.conditional(cond_b, tr_b(u_b), tol) - equation = tr_a(u_a) - K_a * (u_b_padded / K_b) ** n_sorption - deq_du_a = ufl.diff(equation, u_a) * tr_a(v_a) - deq_du_b = ufl.diff(equation, u_b) * tr_b(v_b) - 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_untraced)[0], v_a) + deq_du_b = ufl.inner(ufl.diff(equation, u_b_untraced)[0], v_b) + F_a = ( interface.penalty_term * equation * deq_du_a * dInterface(interface.id) ) - F_1 = ( - -interface.penalty_term + 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_1, n_0) * dInterface( interface.id ) - 0.5 * mixed_term( - v_b, (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_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_1 * 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) ) From aaff7f2e15c7132c4eb34789856ccf2723358055 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 17 Sep 2025 15:08:01 -0400 Subject: [PATCH 3/5] fix subdomain_a --- src/festim/hydrogen_transport_problem.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 2dab95178..d4efd321a 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -1596,6 +1596,8 @@ def mixed_term(u, v, n): 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 From b3746488cc4322be87483d901fdd6a125968ad6f Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 17 Sep 2025 16:30:50 -0400 Subject: [PATCH 4/5] untraced to unrestrcited --- src/festim/hydrogen_transport_problem.py | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index d4efd321a..1b552d5df 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -1551,13 +1551,13 @@ def mixed_term(u, v, n): # use other one for B match subdomain_0.material.solubility_law: case SolubilityLaw.SIEVERT: - u_a, u_a_untraced, v_a, K_a = ( + u_a, u_a_unrestricted, v_a, K_a = ( u_0, u_0_unrestricted, v_0, K_0, ) - u_b, u_b_untraced, v_b, K_b = ( + u_b, u_b_unrestricted, v_b, K_b = ( u_1, u_1_unrestricted, v_1, @@ -1566,13 +1566,13 @@ def mixed_term(u, v, n): subdomain_a = subdomain_0 subdomain_b = subdomain_1 case SolubilityLaw.HENRY: - u_a, u_a_untraced, v_a, K_a = ( + u_a, u_a_unrestricted, v_a, K_a = ( u_1, u_1_unrestricted, v_1, K_1, ) - u_b, u_b_untraced, v_b, K_b = ( + u_b, u_b_unrestricted, v_b, K_b = ( u_0, u_0_unrestricted, v_0, @@ -1583,13 +1583,13 @@ def mixed_term(u, v, n): n_sorption = 0.5 else: - u_a, u_a_untraced, v_a, K_a = ( + u_a, u_a_unrestricted, v_a, K_a = ( u_0, u_0_unrestricted, v_0, K_0, ) - u_b, u_b_untraced, v_b, K_b = ( + u_b, u_b_unrestricted, v_b, K_b = ( u_1, u_1_unrestricted, v_1, @@ -1607,8 +1607,12 @@ def mixed_term(u, v, n): 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_untraced)[0], v_a) - deq_du_b = ufl.inner(ufl.diff(equation, u_b_untraced)[0], v_b) + 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 * equation From b3da603c197f03d3a9a577d4fd9e71c1e7b617cc Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 17 Sep 2025 16:34:10 -0400 Subject: [PATCH 5/5] fixed nitsche --- src/festim/hydrogen_transport_problem.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 1b552d5df..323c432f7 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -1629,10 +1629,10 @@ def mixed_term(u, v, n): subdomain_a.F += F_a subdomain_b.F += F_b case InterfaceMethod.nitsche: - F_0 = -0.5 * mixed_term((u_0 + u_1), v_1, 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_1, (u_0 / K_0 - u_1 / K_1), n_0 + v_0, (u_0 / K_0 - u_1 / K_1), n_0 ) * dInterface(interface.id) F_1 = +0.5 * mixed_term((u_0 + u_1), v_1, n_0) * dInterface( @@ -1645,7 +1645,7 @@ def mixed_term(u, v, n): * gamma / (h_0 + h_1) * (u_0 / K_0 - u_1 / K_1) - * v_1 + * v_0 * dInterface(interface.id) ) F_1 += (