From 853fd0ffb674d1d55d6d92b358eb6ce882805ba2 Mon Sep 17 00:00:00 2001 From: hmonroedd Date: Thu, 9 Apr 2026 17:28:16 -0400 Subject: [PATCH 01/12] added hydrostatic pressure profile --- example2.py | 4 +++- src/sparging/inputs.py | 1 + src/sparging/model.py | 31 ++++++++++++++++++++++++------- test/standard_input.json | 1 + test/test_solve.py | 1 + 5 files changed, 30 insertions(+), 8 deletions(-) diff --git a/example2.py b/example2.py index 3382f5c..1ebead9 100644 --- a/example2.py +++ b/example2.py @@ -21,7 +21,7 @@ geom = ColumnGeometry( area=0.2 * ureg.m**2, - height=1.0 * ureg.m, + height=10 * ureg.m, nozzle_diameter=0.001 * ureg.m, nb_nozzle=10 * ureg.dimensionless, ) @@ -62,6 +62,8 @@ def profile_source_T(z: pint.Quantity): t_final=3 * ureg.days, signal_irr=lambda t: 1 if t < 12 * ureg.hour else 0, signal_sparging=lambda t: 1, + profile_pressure_hydrostatic=False, + # profile_source_T=profile_source_T, ) output = my_simulation.solve() diff --git a/src/sparging/inputs.py b/src/sparging/inputs.py index b61b7e4..efe2b33 100644 --- a/src/sparging/inputs.py +++ b/src/sparging/inputs.py @@ -70,6 +70,7 @@ class SimulationInput: h_l: pint.Quantity K_s: pint.Quantity P_bottom: pint.Quantity + rho_l: pint.Quantity eps_g: pint.Quantity E_g: pint.Quantity D_l: pint.Quantity diff --git a/src/sparging/model.py b/src/sparging/model.py index 94fc5bd..9d1d6b2 100644 --- a/src/sparging/model.py +++ b/src/sparging/model.py @@ -16,7 +16,7 @@ import json from pathlib import Path -from sparging.config import ureg +from sparging.config import ureg, const_g from sparging.inputs import SimulationInput @@ -139,6 +139,13 @@ class Simulation: signal_irr: callable[pint.Quantity] signal_sparging: callable[pint.Quantity] profile_source_T: callable[pint.Quantity] | None = None + profile_pressure_hydrostatic: bool = False + + def hydrostatic_pressure(self, x: pint.Quantity) -> pint.Quantity: + """returns the hydrostatic pressure at a given height x in the tank given P_bottom""" + rho = self.sim_input.rho_l + g = const_g + return (self.sim_input.P_bottom + rho * g * x).to("Pa") def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None): # unpack pint.Quantities @@ -197,24 +204,34 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None else: # homogeneous generation gen_T2 = gen_T2_mag + P = None + if self.profile_pressure_hydrostatic: + P_prof = dolfinx.fem.Function(V_profile) + P_prof.interpolate( + lambda x: self.hydrostatic_pressure(x[0] * ureg.m).magnitude + ) + P = P_prof + else: + P = dolfinx.fem.Constant(mesh, PETSc.ScalarType(P_0)) + # VARIATIONAL FORMULATION # mass transfer rate J_T2 = ( - a * h_l_const * (c_T2 - K_s * (P_0 * y_T2 + EPS)) + a * h_l_const * (c_T2 - K_s * (P * y_T2 + EPS)) ) # TODO pressure shouldn't be a constant (use hydrostatic pressure profile), how to deal with this ? -> use fem.Expression ? F = 0 # variational formulation # transient terms F += eps_l * ((c_T2 - c_T2_n) / dt) * v_c * ufl.dx - F += eps_g * 1 / (const.R * T) * (P_0 * (y_T2 - y_T2_n) / dt) * v_y * ufl.dx + F += eps_g * 1 / (const.R * T) * (P * (y_T2 - y_T2_n) / dt) * v_y * ufl.dx # diffusion/dispersion terms #TODO shouldn't use D_l, transport of T in liquid is dominated by dispersive effects due to gas sparging, find dispersion coeff for steady liquid in gas bubbles F += eps_l * D_l * ufl.dot(ufl.grad(c_T2), ufl.grad(v_c)) * ufl.dx # NOTE remove diffusive term for gas for now for mass balance - # F += eps_g * E_g * ufl.dot(ufl.grad(P_0 * y_T2), ufl.grad(v_y)) * ufl.dx + # F += eps_g * E_g * ufl.dot(ufl.grad(P * y_T2), ufl.grad(v_y)) * ufl.dx # mass exchange (coupling term) F += J_T2 * v_c * ufl.dx - J_T2 * v_y * ufl.dx @@ -226,7 +243,7 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None F += ( 1 / (const.R * T) - * ufl.inner(ufl.dot(ufl.grad(P_0 * y_T2), vel), v_y) + * ufl.inner(ufl.dot(ufl.grad(P * y_T2), vel), v_y) * ufl.dx ) @@ -312,7 +329,7 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None flux_T2 = dolfinx.fem.assemble_scalar( dolfinx.fem.form( - tank_area * vel_x * P_0 / (const.R * T) * y_T2_post * ds(2) + tank_area * vel_x * P / (const.R * T) * y_T2_post * ds(2) ) ) # total T flux at the outlet [mol/s] @@ -330,7 +347,7 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None n = ufl.FacetNormal(mesh) flux_T2_inlet = dolfinx.fem.assemble_scalar( - dolfinx.fem.form(-E_g * ufl.inner(ufl.grad(P_0 * y_T2_post), n) * ds(1)) + dolfinx.fem.form(-E_g * ufl.inner(ufl.grad(P * y_T2_post), n) * ds(1)) ) # total T dispersive flux at the inlet [Pa T2 /s/m2] flux_T2_inlet *= 1 / (const.R * T) * T2_to_T # convert to mol T/s/m2 flux_T2_inlet *= tank_area # convert to mol T/s diff --git a/test/standard_input.json b/test/standard_input.json index 17b74fa..855b9a6 100644 --- a/test/standard_input.json +++ b/test/standard_input.json @@ -7,6 +7,7 @@ "h_l": "2.750e-05 m / s", "K_s": "6.366e-04 mol / m ** 3 / Pa", "P_bottom": "1.208e+05 Pa", + "rho_l": "1.991e+03 kg / m ** 3", "eps_g": "4.091e-04", "E_g": "1.111e-02 m ** 2 / s", "D_l": "2.857e-09 m ** 2 / s", diff --git a/test/test_solve.py b/test/test_solve.py index aceb223..fd3c6f0 100644 --- a/test/test_solve.py +++ b/test/test_solve.py @@ -13,6 +13,7 @@ temperature=600 * ureg.celsius, a=0.5 * ureg("1/m"), h_l=3e-5 * ureg("m/s"), + rho_l=2000 * ureg("kg/m^3"), K_s=1e-4 * ureg("mol/m**3/Pa"), P_bottom=1.2 * ureg.bar, eps_g=0.001 * ureg.dimensionless, From aab8a24b4b0e08320f00e79ae2abca99a0cecb58 Mon Sep 17 00:00:00 2001 From: hmonroedd Date: Mon, 13 Apr 2026 09:25:39 -0400 Subject: [PATCH 02/12] Added J_T2 in SimulationOutput --- example2.py => sparging101.py | 21 +++++++++--------- src/sparging/animation.py | 2 +- src/sparging/model.py | 41 ++++++++++++++++++++++++++--------- 3 files changed, 43 insertions(+), 21 deletions(-) rename example2.py => sparging101.py (80%) diff --git a/example2.py b/sparging101.py similarity index 80% rename from example2.py rename to sparging101.py index 1ebead9..e6cd626 100644 --- a/example2.py +++ b/sparging101.py @@ -21,7 +21,7 @@ geom = ColumnGeometry( area=0.2 * ureg.m**2, - height=10 * ureg.m, + height=1 * ureg.m, nozzle_diameter=0.001 * ureg.m, nb_nozzle=10 * ureg.dimensionless, ) @@ -59,20 +59,21 @@ def profile_source_T(z: pint.Quantity): my_simulation = Simulation( my_input, - t_final=3 * ureg.days, - signal_irr=lambda t: 1 if t < 12 * ureg.hour else 0, + t_final=5 * ureg.days, + signal_irr=lambda t: 1 if t < 8 * ureg.hour else 0, signal_sparging=lambda t: 1, profile_pressure_hydrostatic=False, # profile_source_T=profile_source_T, ) -output = my_simulation.solve() -# # save output to file -# output.profiles_to_csv(f"output_{tank_height}m.csv") +if __name__ == "__main__": + output = my_simulation.solve() -# # plot results -# from sparging import plotting -# plotting.plot_animation(output) + # # save output to file + # output.profiles_to_csv(f"output_{tank_height}m.csv") + # # plot results + # from sparging import plotting + # plotting.plot_animation(output) -animation.create_animation(output, show_activity=True) + animation.create_animation(output, show_activity=False) diff --git a/src/sparging/animation.py b/src/sparging/animation.py index c58ef5a..6baacfd 100644 --- a/src/sparging/animation.py +++ b/src/sparging/animation.py @@ -43,7 +43,7 @@ def __init__( self.x_y = results.x_y self.inventories_T2_salt = results.inventories_T2_salt self.source_T2 = ( - None if results.source_T2 is None else np.array(results.source_T2) + None if results.sources_T2 is None else np.array(results.sources_T2) ) self.fluxes_T2 = ( None if results.fluxes_T2 is None else np.array(results.fluxes_T2) diff --git a/src/sparging/model.py b/src/sparging/model.py index 9d1d6b2..6b66007 100644 --- a/src/sparging/model.py +++ b/src/sparging/model.py @@ -37,20 +37,24 @@ @dataclass -class SimulationResults: +class SimulationResults: # TODO implement pint in this class # TODO change list to np.array times: list c_T2_solutions: list y_T2_solutions: list + J_T2_solutions: list x_ct: np.ndarray x_y: np.ndarray inventories_T2_salt: np.ndarray - source_T2: list + sources_T2: list fluxes_T2: list sim_input: SimulationInput + dt: int | None = None + dx: int | None = None keys_to_ignore_results = [ # TODO do it the other way: keys_to_include_results "c_T2_solutions", "y_T2_solutions", + "J_T2_solutions", "x_ct", "x_y", "inventories_T2_salt", @@ -200,19 +204,28 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None gen_T2_prof.interpolate( lambda x: x[0] * 0 + self.profile_source_T(x[0] * ureg.m) ) - gen_T2 = gen_T2_mag * gen_T2_prof + integral_gen_T2 = dolfinx.fem.assemble_scalar( + dolfinx.fem.form(gen_T2_prof * ufl.dx) + ) else: # homogeneous generation - gen_T2 = gen_T2_mag + gen_T2_prof = dolfinx.fem.Constant( + mesh, PETSc.ScalarType(1.0 / tank_height) + ) + integral_gen_T2 = dolfinx.fem.Constant(mesh, PETSc.ScalarType(1.0)) + + gen_T2_prof_normalized = gen_T2_prof / integral_gen_T2 - P = None + gen_T2 = gen_T2_mag * gen_T2_prof_normalized + + P_prof = dolfinx.fem.Function(V_profile) if self.profile_pressure_hydrostatic: - P_prof = dolfinx.fem.Function(V_profile) P_prof.interpolate( lambda x: self.hydrostatic_pressure(x[0] * ureg.m).magnitude ) - P = P_prof else: - P = dolfinx.fem.Constant(mesh, PETSc.ScalarType(P_0)) + P_prof.interpolate(lambda x: x[0] * 0 + P_0) + + P = P_prof # VARIATIONAL FORMULATION @@ -296,6 +309,7 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None times = [] c_T2_solutions = [] y_T2_solutions = [] + J_T2_solutions = [] sources_T2 = [] fluxes_T2 = [] inventories_T2_salt = [] @@ -303,9 +317,9 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None # SOLVE t = 0 while t < t_final: + # update time-dependent terms gen_T2_mag.value = source_T2 * self.signal_irr(t * ureg.s) h_l_const.value = h_l * self.signal_sparging(t * ureg.s) - """ utiliser ufl.conditional TODO""" problem.solve() @@ -317,12 +331,16 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None c_T2_vals = u.x.array[ct_dofs][ct_sort_coords] y_T2_vals = u.x.array[y_dofs][y_sort_coords] + J_T2_vals = ( + a * h_l_const.value * (c_T2_vals - K_s * (P.x.array * y_T2_vals + EPS)) + ) # TODO there is a clever way of getting J_T2 for sure # store time and solution # TODO give units to the results times.append(t) c_T2_solutions.append(c_T2_vals.copy()) y_T2_solutions.append(y_T2_vals.copy()) + J_T2_solutions.append(J_T2_vals.copy()) sources_T2.append( source_T2 * self.signal_irr(t * ureg.s) * tank_volume ) # total T generation rate in the tank [mol/s] TODO useless: signal_irr is already given @@ -372,11 +390,14 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None times=times, c_T2_solutions=c_T2_solutions, y_T2_solutions=y_T2_solutions, + J_T2_solutions=J_T2_solutions, x_ct=x_ct, x_y=x_y, inventories_T2_salt=inventories_T2_salt, - source_T2=sources_T2, + sources_T2=sources_T2, fluxes_T2=fluxes_T2, sim_input=self.sim_input, + dt=dt, + dx=dx, ) return results From a7f442c8bf0f8976bdfacede5877215f799b5eab Mon Sep 17 00:00:00 2001 From: hmonroedd Date: Mon, 13 Apr 2026 11:26:24 -0400 Subject: [PATCH 03/12] Added liquid phase dispersion E_l + updated test -> got a dispersion problem, need to fix BCs --- sparging101.py | 4 ++-- src/sparging/correlations.py | 25 ++++++++++++++++--------- src/sparging/inputs.py | 2 ++ src/sparging/model.py | 9 ++++----- test/standard_input.json | 1 + test/test_solve.py | 1 + 6 files changed, 26 insertions(+), 16 deletions(-) diff --git a/sparging101.py b/sparging101.py index e6cd626..8d1e8b3 100644 --- a/sparging101.py +++ b/sparging101.py @@ -60,7 +60,7 @@ def profile_source_T(z: pint.Quantity): my_simulation = Simulation( my_input, t_final=5 * ureg.days, - signal_irr=lambda t: 1 if t < 8 * ureg.hour else 0, + signal_irr=lambda t: 1 if t < 3 * ureg.days else 0, signal_sparging=lambda t: 1, profile_pressure_hydrostatic=False, # profile_source_T=profile_source_T, @@ -76,4 +76,4 @@ def profile_source_T(z: pint.Quantity): # from sparging import plotting # plotting.plot_animation(output) - animation.create_animation(output, show_activity=False) + animation.create_animation(output, show_activity=True) diff --git a/src/sparging/correlations.py b/src/sparging/correlations.py index 9f8f09b..a586593 100644 --- a/src/sparging/correlations.py +++ b/src/sparging/correlations.py @@ -30,6 +30,7 @@ class CorrelationType(enum.Enum): # TODO do we really use it ? REYNOLDS_NUMBER = "Re" BUBBLE_VELOCITY = "u_g0" GAS_PHASE_DISPERSION = "E_g" + LIQUID_PHASE_DISPERSION = "E_l" PRESSURE = "P" FLOW_RATE = "flow_g_mol" INTERFACIAL_AREA = "a" @@ -306,10 +307,23 @@ def __contains__(self, key: str | Correlation): ) all_correlations.append(h_l_briggs) +E_l = Correlation( + identifier="E_l", + function=lambda tank_diameter, u_g0: ureg.Quantity( + 0.678 * tank_diameter.magnitude**1.4 * u_g0.magnitude**0.3, "m**2/s" + ), # liquid phase axial dispersion coefficient + corr_type=CorrelationType.LIQUID_PHASE_DISPERSION, + source="Deckwer 1974", + description="liquid phase axial dispersion coefficient, assumed equal to diffusivity of tritium in liquid FLiBe", + input_units=["m", "m/s"], + output_units="m**2/s", +) +all_correlations.append(E_l) + E_g = Correlation( identifier="E_g", - function=lambda tank_diameter, u_g0: get_E_g( - diameter=tank_diameter, u_g=u_g0 + function=lambda tank_diameter, u_g0: ( + 0.2 * ureg("1/m") * tank_diameter**2 * u_g0 ), # gas phase axial dispersion coefficient corr_type=CorrelationType.GAS_PHASE_DISPERSION, source="Malara 1995", @@ -474,10 +488,3 @@ def get_h_briggs(Re: float, Sc: float, D_l: float, d_b: float) -> float: Sh = 0.089 * Re**0.69 * Sc**0.33 # Sherwood number h_l = Sh * D_l / d_b return h_l - - -def get_E_g(diameter: float, u_g: float) -> float: - """gas phase axial dispersion coefficient [m2/s], Malara 1995 correlation - models dispersion of the gas velocity distribution around the mean bubble velocity""" - E_g = 0.2 * ureg("1/m") * diameter**2 * u_g - return E_g diff --git a/src/sparging/inputs.py b/src/sparging/inputs.py index efe2b33..8ac91dc 100644 --- a/src/sparging/inputs.py +++ b/src/sparging/inputs.py @@ -57,6 +57,7 @@ class SpargingParameters: d_b: pint.Quantity | Correlation | None = None rho_g: pint.Quantity | Correlation | None = None E_g: pint.Quantity | Correlation | None = None + E_l: pint.Quantity | Correlation | None = None a: pint.Quantity | Correlation | None = None @@ -73,6 +74,7 @@ class SimulationInput: rho_l: pint.Quantity eps_g: pint.Quantity E_g: pint.Quantity + E_l: pint.Quantity D_l: pint.Quantity source_T: pint.Quantity diff --git a/src/sparging/model.py b/src/sparging/model.py index 6b66007..484c118 100644 --- a/src/sparging/model.py +++ b/src/sparging/model.py @@ -164,7 +164,8 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None T = self.sim_input.temperature.to("K").magnitude eps_g = self.sim_input.eps_g.to("dimensionless").magnitude E_g = self.sim_input.E_g.to("m**2/s").magnitude - D_l = self.sim_input.D_l.to("m**2/s").magnitude + E_l = self.sim_input.E_l.to("m**2/s").magnitude + D_l = self.sim_input.D_l.to("m**2/s").magnitude # not needed (included in h_l) u_g0 = self.sim_input.u_g0.to("m/s").magnitude source_T2 = self.sim_input.source_T.to("molT2/s/m**3").magnitude @@ -240,10 +241,8 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None F += eps_l * ((c_T2 - c_T2_n) / dt) * v_c * ufl.dx F += eps_g * 1 / (const.R * T) * (P * (y_T2 - y_T2_n) / dt) * v_y * ufl.dx - # diffusion/dispersion terms #TODO shouldn't use D_l, transport of T in liquid is dominated by dispersive effects due to gas sparging, find dispersion coeff for steady liquid in gas bubbles - F += eps_l * D_l * ufl.dot(ufl.grad(c_T2), ufl.grad(v_c)) * ufl.dx - - # NOTE remove diffusive term for gas for now for mass balance + # dispersive terms NOTE removed dispersive terms for now for mass balance + # F += eps_l * E_l * ufl.dot(ufl.grad(c_T2), ufl.grad(v_c)) * ufl.dx # F += eps_g * E_g * ufl.dot(ufl.grad(P * y_T2), ufl.grad(v_y)) * ufl.dx # mass exchange (coupling term) diff --git a/test/standard_input.json b/test/standard_input.json index 855b9a6..77f4c5f 100644 --- a/test/standard_input.json +++ b/test/standard_input.json @@ -10,6 +10,7 @@ "rho_l": "1.991e+03 kg / m ** 3", "eps_g": "4.091e-04", "E_g": "1.111e-02 m ** 2 / s", + "E_l": "1.648e-01 m ** 2 / s", "D_l": "2.857e-09 m ** 2 / s", "source_T": "8.303e-16 molT / m ** 3 / s" } \ No newline at end of file diff --git a/test/test_solve.py b/test/test_solve.py index fd3c6f0..2894840 100644 --- a/test/test_solve.py +++ b/test/test_solve.py @@ -18,6 +18,7 @@ P_bottom=1.2 * ureg.bar, eps_g=0.001 * ureg.dimensionless, E_g=1e-2 * ureg("m^2/s"), + E_l=1e-1 * ureg("m^2/s"), D_l=3e-9 * ureg("m^2/s"), source_T=8e-16 * ureg("molT/m^3/s"), ) From 42ea13b5f30a09d3b3ea50475874fdf235bbc79e Mon Sep 17 00:00:00 2001 From: hmonroedd Date: Mon, 13 Apr 2026 16:38:24 -0400 Subject: [PATCH 04/12] fixed inventory dependance on tank height (missing factor when integrating) --- src/sparging/model.py | 35 ++++++++++++++++++++++------------- 1 file changed, 22 insertions(+), 13 deletions(-) diff --git a/src/sparging/model.py b/src/sparging/model.py index 484c118..7ad8bec 100644 --- a/src/sparging/model.py +++ b/src/sparging/model.py @@ -216,7 +216,9 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None gen_T2_prof_normalized = gen_T2_prof / integral_gen_T2 - gen_T2 = gen_T2_mag * gen_T2_prof_normalized + gen_T2 = ( + gen_T2_mag * gen_T2_prof_normalized + ) # so that integral over volume of gen_T2 = source_T2 P_prof = dolfinx.fem.Function(V_profile) if self.profile_pressure_hydrostatic: @@ -270,12 +272,12 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None dolfinx.fem.Constant(mesh, 0.0), dolfinx.fem.locate_dofs_topological(V.sub(1), fdim, gas_inlet_facets), V.sub(1), - ) + ) # y_T2 = 0 at gas inlet bc2 = dolfinx.fem.dirichletbc( dolfinx.fem.Constant(mesh, 0.0), dolfinx.fem.locate_dofs_topological(V.sub(0), fdim, gas_outlet_facets), V.sub(0), - ) + ) # c_T2 = 0 at gas outlet TODO should apply Neumann BC rather # Custom measure all_facets = np.concatenate((gas_inlet_facets, gas_outlet_facets)) @@ -344,12 +346,18 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None source_T2 * self.signal_irr(t * ureg.s) * tank_volume ) # total T generation rate in the tank [mol/s] TODO useless: signal_irr is already given + n = ufl.FacetNormal(mesh) + + # flux_T2 = dolfinx.fem.assemble_scalar( + # dolfinx.fem.form( + # eps_g * vel_x * P / (const.R * T) * y_T2_post * tank_area * ds(2) + # ) + # ) # total T flux at the outlet [mol/s] flux_T2 = dolfinx.fem.assemble_scalar( dolfinx.fem.form( - tank_area * vel_x * P / (const.R * T) * y_T2_post * ds(2) + vel_x * P / (const.R * T) * y_T2_post * tank_area * ds(2) ) - ) # total T flux at the outlet [mol/s] - + ) # flux_T_inlet = dolfinx.fem.assemble_scalar( # dolfinx.fem.form( # tank_area @@ -362,18 +370,19 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None # ) # ) # total T dispersive flux at the inlet [mol/s] - n = ufl.FacetNormal(mesh) - flux_T2_inlet = dolfinx.fem.assemble_scalar( - dolfinx.fem.form(-E_g * ufl.inner(ufl.grad(P * y_T2_post), n) * ds(1)) - ) # total T dispersive flux at the inlet [Pa T2 /s/m2] - flux_T2_inlet *= 1 / (const.R * T) * T2_to_T # convert to mol T/s/m2 - flux_T2_inlet *= tank_area # convert to mol T/s + # flux_T2_inlet = dolfinx.fem.assemble_scalar( + # dolfinx.fem.form(-E_g * ufl.inner(ufl.grad(P * y_T2_post), n) * ds(1)) + # ) # total T dispersive flux at the inlet [Pa T2 /s/m2] + # flux_T2_inlet *= 1 / (const.R * T) * T2_to_T # convert to mol T/s/m2 + # flux_T2_inlet *= tank_area # convert to mol T/s # fluxes_T2.append(flux_T2 + flux_T2_inlet) fluxes_T2.append(flux_T2) inventory_T2_salt = dolfinx.fem.assemble_scalar( - dolfinx.fem.form(c_T2_post * ufl.dx) + dolfinx.fem.form( + c_T2_post * tank_height * ufl.dx + ) # NOTE looks like dx € [0,1], inventory should be constant no matter the height when specifying tbr and n_gen_rate ) inventory_T2_salt *= tank_area # get total amount of T2 in [mol] inventories_T2_salt.append(inventory_T2_salt) From d6d61300e86323ee7125c3736ced5d64ef91fd11 Mon Sep 17 00:00:00 2001 From: hmonroedd Date: Mon, 13 Apr 2026 17:40:58 -0400 Subject: [PATCH 05/12] fixed inventory dependance on tank height (truly this time) by making distinction between normalized and integral source_T --- sparging101.py | 22 ++++++++++++++-------- src/sparging/correlations.py | 6 +++--- src/sparging/inputs.py | 8 ++++++-- src/sparging/model.py | 31 +++++++++++-------------------- 4 files changed, 34 insertions(+), 33 deletions(-) diff --git a/sparging101.py b/sparging101.py index 8d1e8b3..a9a8452 100644 --- a/sparging101.py +++ b/sparging101.py @@ -21,7 +21,7 @@ geom = ColumnGeometry( area=0.2 * ureg.m**2, - height=1 * ureg.m, + height=0.1 * ureg.m, nozzle_diameter=0.001 * ureg.m, nb_nozzle=10 * ureg.dimensionless, ) @@ -42,26 +42,28 @@ h_l=all_correlations("h_l_briggs"), ) - # class method from_parameters that takes in objects like ColumnGeometry, BreederMaterial, OperatingParameters and returns a SimulationInput object with the appropriate correlations for the given parameters. This method should be able to handle cases where some of the parameters are already provided as correlations and should not overwrite them. my_input = SimulationInput.from_parameters( geom, flibe, operating_params, sparging_params ) logger.info(my_input) +print(my_input.source_T_norm) +print(my_input.source_T_int.to("molT/s")) + def profile_source_T(z: pint.Quantity): import numpy as np - # return np.sin(np.pi / (1 * ureg.m) * z) - return 0.5 * (1 + np.cos(0.5 * np.pi / (1 * ureg.m) * z)) + return np.sin(np.pi / (1 * ureg.m) * z) + # return 0.5 * (1 + np.cos(0.5 * np.pi / (1 * ureg.m) * z)) my_simulation = Simulation( my_input, - t_final=5 * ureg.days, - signal_irr=lambda t: 1 if t < 3 * ureg.days else 0, - signal_sparging=lambda t: 1, + t_final=200 * ureg.s, + signal_irr=lambda t: 1 if t < 100 * ureg.s else 0, + signal_sparging=lambda t: 0, profile_pressure_hydrostatic=False, # profile_source_T=profile_source_T, ) @@ -76,4 +78,8 @@ def profile_source_T(z: pint.Quantity): # from sparging import plotting # plotting.plot_animation(output) - animation.create_animation(output, show_activity=True) + # import matplotlib.pyplot as plt + + # plt.plot(output.times, output.inventories_T2_salt) + # plt.show() + # animation.create_animation(output, show_activity=False) diff --git a/src/sparging/correlations.py b/src/sparging/correlations.py index a586593..6dd8c59 100644 --- a/src/sparging/correlations.py +++ b/src/sparging/correlations.py @@ -393,8 +393,8 @@ def __contains__(self, key: str | Correlation): ) all_correlations.append(specific_interfacial_area) -source_T_from_tbr = Correlation( - identifier="source_T", +source_T_normalized = Correlation( + identifier="source_T_norm", function=lambda tbr, n_gen_rate, tank_volume: ( tbr * n_gen_rate / tank_volume ), # source term for tritium generation calculated from TBR and neutron generation rate @@ -402,7 +402,7 @@ def __contains__(self, key: str | Correlation): input_units=["triton/neutron", "neutron/s", "m**3"], output_units="molT/m**3/s", ) -all_correlations.append(source_T_from_tbr) +all_correlations.append(source_T_normalized) def get_d_b( diff --git a/src/sparging/inputs.py b/src/sparging/inputs.py index 8ac91dc..a545237 100644 --- a/src/sparging/inputs.py +++ b/src/sparging/inputs.py @@ -46,7 +46,7 @@ class OperatingParameters: P_bottom: pint.Quantity | Correlation | None = None tbr: pint.Quantity | None = None n_gen_rate: pint.Quantity | None = None - source_T: pint.Quantity | Correlation | None = None + source_T_norm: pint.Quantity | Correlation | None = None @dataclass @@ -76,12 +76,16 @@ class SimulationInput: E_g: pint.Quantity E_l: pint.Quantity D_l: pint.Quantity - source_T: pint.Quantity + source_T_norm: pint.Quantity @property def volume(self): return self.area * self.height + @property + def source_T_int(self): + return self.source_T_norm * self.volume + def __post_init__(self): # make sure there are only pint.Quantity or callables in the input, otherwise raise an error for key in self.__dataclass_fields__.keys(): diff --git a/src/sparging/model.py b/src/sparging/model.py index 7ad8bec..82e9c0b 100644 --- a/src/sparging/model.py +++ b/src/sparging/model.py @@ -167,7 +167,7 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None E_l = self.sim_input.E_l.to("m**2/s").magnitude D_l = self.sim_input.D_l.to("m**2/s").magnitude # not needed (included in h_l) u_g0 = self.sim_input.u_g0.to("m/s").magnitude - source_T2 = self.sim_input.source_T.to("molT2/s/m**3").magnitude + source_T2_norm = self.sim_input.source_T_norm.to("molT2/s/m**3").magnitude dt = dt.to("seconds").magnitude if dt is not None else t_final / 1000 dx = dx.to("m").magnitude if dx is not None else tank_height / 1000 @@ -196,29 +196,22 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None h_l_const = dolfinx.fem.Constant(mesh, PETSc.ScalarType(h_l)) - gen_T2_mag = dolfinx.fem.Constant( - mesh, source_T2 * self.signal_irr(0 * ureg.s) + gen_T2_norm = dolfinx.fem.Constant( + mesh, source_T2_norm * self.signal_irr(0 * ureg.s) ) # magnitude of the generation term - gen_T2 = None + if self.profile_source_T is not None: # spatially varying profile is provided gen_T2_prof = dolfinx.fem.Function(V_profile) gen_T2_prof.interpolate( lambda x: x[0] * 0 + self.profile_source_T(x[0] * ureg.m) ) - integral_gen_T2 = dolfinx.fem.assemble_scalar( - dolfinx.fem.form(gen_T2_prof * ufl.dx) - ) else: # homogeneous generation - gen_T2_prof = dolfinx.fem.Constant( - mesh, PETSc.ScalarType(1.0 / tank_height) - ) - integral_gen_T2 = dolfinx.fem.Constant(mesh, PETSc.ScalarType(1.0)) - - gen_T2_prof_normalized = gen_T2_prof / integral_gen_T2 + gen_T2_prof = dolfinx.fem.Constant(mesh, PETSc.ScalarType(1.0)) gen_T2 = ( - gen_T2_mag * gen_T2_prof_normalized - ) # so that integral over volume of gen_T2 = source_T2 + gen_T2_norm * gen_T2_prof + ) # gen_T2_norm is already mormalized by tank volume + # user can choose to specify a normalized profile_source_T so that mean(gen_T2) = gen_T2_norm = 0.5 * source_T_norm P_prof = dolfinx.fem.Function(V_profile) if self.profile_pressure_hydrostatic: @@ -319,7 +312,7 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None t = 0 while t < t_final: # update time-dependent terms - gen_T2_mag.value = source_T2 * self.signal_irr(t * ureg.s) + gen_T2_norm.value = source_T2_norm * self.signal_irr(t * ureg.s) h_l_const.value = h_l * self.signal_sparging(t * ureg.s) problem.solve() @@ -343,7 +336,7 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None y_T2_solutions.append(y_T2_vals.copy()) J_T2_solutions.append(J_T2_vals.copy()) sources_T2.append( - source_T2 * self.signal_irr(t * ureg.s) * tank_volume + source_T2_norm * self.signal_irr(t * ureg.s) * tank_volume ) # total T generation rate in the tank [mol/s] TODO useless: signal_irr is already given n = ufl.FacetNormal(mesh) @@ -380,9 +373,7 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None fluxes_T2.append(flux_T2) inventory_T2_salt = dolfinx.fem.assemble_scalar( - dolfinx.fem.form( - c_T2_post * tank_height * ufl.dx - ) # NOTE looks like dx € [0,1], inventory should be constant no matter the height when specifying tbr and n_gen_rate + dolfinx.fem.form(c_T2_post * ufl.dx) ) inventory_T2_salt *= tank_area # get total amount of T2 in [mol] inventories_T2_salt.append(inventory_T2_salt) From 5c45de6091e1322a7e40878c24aba1bd61970d72 Mon Sep 17 00:00:00 2001 From: hmonroedd Date: Mon, 13 Apr 2026 17:50:20 -0400 Subject: [PATCH 06/12] LIBRA 1L model --- sparging_LIBRA1L.py | 76 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 76 insertions(+) create mode 100644 sparging_LIBRA1L.py diff --git a/sparging_LIBRA1L.py b/sparging_LIBRA1L.py new file mode 100644 index 0000000..0fa2c77 --- /dev/null +++ b/sparging_LIBRA1L.py @@ -0,0 +1,76 @@ +from sparging.config import ureg +from sparging import all_correlations +from sparging import animation +from sparging.model import Simulation +from sparging.inputs import ( + ColumnGeometry, + BreederMaterial, + OperatingParameters, + SpargingParameters, + SimulationInput, +) +import logging +from typing import TYPE_CHECKING + +if TYPE_CHECKING: + import pint + +logger = logging.getLogger(__name__) +logging.basicConfig(level=logging.WARNING) + +geom = ColumnGeometry( + area=170 * ureg.cm**2, + height=7 * ureg.cm, + nozzle_diameter=1.4 * ureg.mm, + nb_nozzle=1 * ureg.dimensionless, +) + +flibe = BreederMaterial( + name="FLiBe", +) + +operating_params = OperatingParameters( + temperature=600 * ureg.celsius, + P_top=1 * ureg.atm, + flow_g_mol=40 * ureg.sccm, + tbr=2e-3 * ureg("triton / neutron"), + n_gen_rate=1e9 * ureg("neutron / s"), +) + +sparging_params = SpargingParameters( + h_l=all_correlations("h_l_briggs"), +) + +# class method from_parameters that takes in objects like ColumnGeometry, BreederMaterial, OperatingParameters and returns a SimulationInput object with the appropriate correlations for the given parameters. This method should be able to handle cases where some of the parameters are already provided as correlations and should not overwrite them. +my_input = SimulationInput.from_parameters( + geom, flibe, operating_params, sparging_params +) +logger.info(my_input) +print(my_input.source_T_norm) +print(my_input.source_T_int.to("molT/s")) + +n_fluence = 2.5e13 * ureg("neutron") +n_gen_rate = operating_params.n_gen_rate +t_irr = n_fluence / n_gen_rate +print(f"t_irr = {t_irr.to('seconds')}") +my_simulation = Simulation( + my_input, + t_final=50 * ureg.days, + signal_irr=lambda t: 1 if t < t_irr else 0, + signal_sparging=lambda t: 0 if t < t_irr else 1, + # signal_sparging=lambda t: 0, + profile_pressure_hydrostatic=False, + profile_source_T=lambda z: 1, +) + +if __name__ == "__main__": + output = my_simulation.solve() + + # # save output to file + # output.profiles_to_csv(f"output_{tank_height}m.csv") + + # # plot results + # from sparging import plotting + # plotting.plot_animation(output) + + animation.create_animation(output, show_activity=False) From 901c0950f16620816fa8c3e1959f78b6c5341b5e Mon Sep 17 00:00:00 2001 From: hmonroedd Date: Wed, 15 Apr 2026 15:07:27 -0400 Subject: [PATCH 07/12] re-introduced dispersion terms, with neumann BC, mass conservation should be satisfied now --- src/sparging/model.py | 48 +++++++++++++++++++++++++++---------------- 1 file changed, 30 insertions(+), 18 deletions(-) diff --git a/src/sparging/model.py b/src/sparging/model.py index 82e9c0b..89bfdbd 100644 --- a/src/sparging/model.py +++ b/src/sparging/model.py @@ -203,11 +203,12 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None if self.profile_source_T is not None: # spatially varying profile is provided gen_T2_prof = dolfinx.fem.Function(V_profile) gen_T2_prof.interpolate( - lambda x: x[0] * 0 + self.profile_source_T(x[0] * ureg.m) + lambda x: x[0] * 0 + self.profile_source_T(x[0] / tank_height) ) else: # homogeneous generation gen_T2_prof = dolfinx.fem.Constant(mesh, PETSc.ScalarType(1.0)) + # if self.sim_input.normalize_source_T: gen_T2 = ( gen_T2_norm * gen_T2_prof ) # gen_T2_norm is already mormalized by tank volume @@ -236,9 +237,17 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None F += eps_l * ((c_T2 - c_T2_n) / dt) * v_c * ufl.dx F += eps_g * 1 / (const.R * T) * (P * (y_T2 - y_T2_n) / dt) * v_y * ufl.dx - # dispersive terms NOTE removed dispersive terms for now for mass balance - # F += eps_l * E_l * ufl.dot(ufl.grad(c_T2), ufl.grad(v_c)) * ufl.dx - # F += eps_g * E_g * ufl.dot(ufl.grad(P * y_T2), ufl.grad(v_y)) * ufl.dx + # dispersive terms + F += eps_l * E_l * ufl.dot(ufl.grad(c_T2), ufl.grad(v_c)) * ufl.dx + F += ( + eps_g + * E_g + * 1 + / (const.R * T) + * ufl.dot(ufl.grad(P * y_T2), ufl.grad(v_y)) + * ufl.dx + ) + # there was missing 1/RT # mass exchange (coupling term) F += J_T2 * v_c * ufl.dx - J_T2 * v_y * ufl.dx @@ -266,11 +275,11 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None dolfinx.fem.locate_dofs_topological(V.sub(1), fdim, gas_inlet_facets), V.sub(1), ) # y_T2 = 0 at gas inlet - bc2 = dolfinx.fem.dirichletbc( - dolfinx.fem.Constant(mesh, 0.0), - dolfinx.fem.locate_dofs_topological(V.sub(0), fdim, gas_outlet_facets), - V.sub(0), - ) # c_T2 = 0 at gas outlet TODO should apply Neumann BC rather + # bc2 = dolfinx.fem.dirichletbc( + # dolfinx.fem.Constant(mesh, 0.0), + # dolfinx.fem.locate_dofs_topological(V.sub(0), fdim, gas_outlet_facets), + # V.sub(0), + # ) # c_T2 = 0 at gas outlet # Custom measure all_facets = np.concatenate((gas_inlet_facets, gas_outlet_facets)) @@ -284,7 +293,8 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None problem = NonlinearProblem( F, u, - bcs=[bc1, bc2], + # bcs=[bc1, bc2], + bcs=[bc1], # Neumann BCs on c_T2 at inlet and outlet are naturally enforced petsc_options_prefix="librasparge", # petsc_options={"snes_monitor": None}, ) @@ -350,7 +360,7 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None dolfinx.fem.form( vel_x * P / (const.R * T) * y_T2_post * tank_area * ds(2) ) - ) + ) # TODO replace with integral of J over volume # flux_T_inlet = dolfinx.fem.assemble_scalar( # dolfinx.fem.form( # tank_area @@ -363,14 +373,16 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None # ) # ) # total T dispersive flux at the inlet [mol/s] - # flux_T2_inlet = dolfinx.fem.assemble_scalar( - # dolfinx.fem.form(-E_g * ufl.inner(ufl.grad(P * y_T2_post), n) * ds(1)) - # ) # total T dispersive flux at the inlet [Pa T2 /s/m2] - # flux_T2_inlet *= 1 / (const.R * T) * T2_to_T # convert to mol T/s/m2 - # flux_T2_inlet *= tank_area # convert to mol T/s + flux_T2_inlet = dolfinx.fem.assemble_scalar( + dolfinx.fem.form( + -eps_g * E_g * ufl.inner(ufl.grad(P * y_T2_post), n) * ds(1) + ) + ) # total T dispersive flux at the inlet [Pa T2 /s/m2] + flux_T2_inlet *= 1 / (const.R * T) # mol T2/s/m2 + flux_T2_inlet *= tank_area # convert to molT2/s - # fluxes_T2.append(flux_T2 + flux_T2_inlet) - fluxes_T2.append(flux_T2) + fluxes_T2.append(flux_T2 + flux_T2_inlet) + # fluxes_T2.append(flux_T2) inventory_T2_salt = dolfinx.fem.assemble_scalar( dolfinx.fem.form(c_T2_post * ufl.dx) From 6bdae661a1cd82f10db53f61ed4782feefc0cb45 Mon Sep 17 00:00:00 2001 From: hmonroedd Date: Wed, 15 Apr 2026 17:05:51 -0400 Subject: [PATCH 08/12] added quantities of the analytical model in SimulationInput --- sparging101.py | 32 +++++++++++++++++++++++++------- src/sparging/correlations.py | 2 +- src/sparging/inputs.py | 31 ++++++++++++++++++++++++++++++- src/sparging/model.py | 1 - 4 files changed, 56 insertions(+), 10 deletions(-) diff --git a/sparging101.py b/sparging101.py index a9a8452..2f13823 100644 --- a/sparging101.py +++ b/sparging101.py @@ -11,6 +11,7 @@ ) import logging from typing import TYPE_CHECKING +import numpy as np if TYPE_CHECKING: import pint @@ -21,7 +22,7 @@ geom = ColumnGeometry( area=0.2 * ureg.m**2, - height=0.1 * ureg.m, + height=2 * ureg.m, nozzle_diameter=0.001 * ureg.m, nb_nozzle=10 * ureg.dimensionless, ) @@ -50,25 +51,42 @@ print(my_input.source_T_norm) print(my_input.source_T_int.to("molT/s")) +print( + f"Concentration at steady state (no dispersion, no PP limited): { + my_input.get_c_T2_SS().to('molT2/m^3') + }" +) +T_99 = (np.log(100) * my_input.get_tau()).to("seconds") +print(f"T_99% = {T_99.to('hours')}") +print(f"Partial pressure number PP = {my_input.get_PP_number()}") -def profile_source_T(z: pint.Quantity): +def profile_source_T(z: pint.Quantity | list[float], height: pint.Quantity = None): import numpy as np - return np.sin(np.pi / (1 * ureg.m) * z) + if isinstance(z, (float, np.ndarray, list)): # non-dimensional height (0 to 1) + return np.pi / 2 * np.sin(np.pi * z) # normalized + if isinstance(z, ureg.Quantity): + if height is None: + raise ValueError("Must provide height if z is a dimensional quantity") + return np.pi / 2 * np.sin(np.pi / height * z) # normalized + else: + raise NotImplementedError("z must be either a float or a pint.Quantity") # return 0.5 * (1 + np.cos(0.5 * np.pi / (1 * ureg.m) * z)) my_simulation = Simulation( my_input, - t_final=200 * ureg.s, - signal_irr=lambda t: 1 if t < 100 * ureg.s else 0, - signal_sparging=lambda t: 0, + t_final=2 * T_99, + signal_irr=lambda t: 1 if t < T_99 else 0, + signal_sparging=lambda t: 1, profile_pressure_hydrostatic=False, # profile_source_T=profile_source_T, ) if __name__ == "__main__": + # my_simulation.sim_input.E_g *= 1e5 + # my_simulation.sim_input.E_l *= 1e-5 output = my_simulation.solve() # # save output to file @@ -82,4 +100,4 @@ def profile_source_T(z: pint.Quantity): # plt.plot(output.times, output.inventories_T2_salt) # plt.show() - # animation.create_animation(output, show_activity=False) + animation.create_animation(output, show_activity=False) diff --git a/src/sparging/correlations.py b/src/sparging/correlations.py index 6dd8c59..7d41f32 100644 --- a/src/sparging/correlations.py +++ b/src/sparging/correlations.py @@ -393,7 +393,7 @@ def __contains__(self, key: str | Correlation): ) all_correlations.append(specific_interfacial_area) -source_T_normalized = Correlation( +source_T_normalized = Correlation( # TODO remettre source_T identifier="source_T_norm", function=lambda tbr, n_gen_rate, tank_volume: ( tbr * n_gen_rate / tank_volume diff --git a/src/sparging/inputs.py b/src/sparging/inputs.py index a545237..29efe89 100644 --- a/src/sparging/inputs.py +++ b/src/sparging/inputs.py @@ -5,6 +5,7 @@ import inspect import numpy as np import logging +from sparging.config import const_R logger = logging.getLogger(__name__) @@ -77,15 +78,43 @@ class SimulationInput: E_l: pint.Quantity D_l: pint.Quantity source_T_norm: pint.Quantity + # normalize_source_T: bool = True @property def volume(self): return self.area * self.height @property - def source_T_int(self): + def eps_l(self): + return 1 - self.eps_g + + @property # TODO enlever, inutile + def source_T_int(self): # rename Q_T2 return self.source_T_norm * self.volume + def set_S_T2(self, val: pint.Quantity): # TODO to use + self.Q_T2 = (val.to("molT2/m**3/s") * self.volume).to("molT2/s") + + def get_tau(self) -> pint.Quantity: + """characteristic time of the sparger under the small partial pressure (SPP) approximation""" + return (self.eps_l / (self.h_l * self.a)).to("seconds") + + def get_c_T2_SS(self) -> pint.Quantity: + return (self.source_T_int / self.volume * 1 / (self.h_l * self.a)).to( + "molT2/m^3" + ) + + def get_PP_number(self) -> pint.Quantity: + """Partial pressure number, ratio of the equivalent T concentration at liquid boundary to the bulk liquid concentration. + If PP << 1, then we are in the small partial pressure (SPP) regime, and if PP ~= 1, then we are in the partial pressure limited (PPL) regime. + """ + return ( + self.K_s + * (const_R * self.temperature) + * self.height + / (self.get_tau() * self.u_g0) + ).to("dimensionless") + def __post_init__(self): # make sure there are only pint.Quantity or callables in the input, otherwise raise an error for key in self.__dataclass_fields__.keys(): diff --git a/src/sparging/model.py b/src/sparging/model.py index 89bfdbd..8ae15cc 100644 --- a/src/sparging/model.py +++ b/src/sparging/model.py @@ -247,7 +247,6 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None * ufl.dot(ufl.grad(P * y_T2), ufl.grad(v_y)) * ufl.dx ) - # there was missing 1/RT # mass exchange (coupling term) F += J_T2 * v_c * ufl.dx - J_T2 * v_y * ufl.dx From 28504606c048be8e7c06710e5ae5cacd1775a943 Mon Sep 17 00:00:00 2001 From: hmonroedd Date: Thu, 16 Apr 2026 11:30:24 -0400 Subject: [PATCH 09/12] source normalization now based on integral source + added associated test + fast_run and dispersion_run options in model.py --- sparging101.py | 19 ++++++---- sparging_LIBRA1L.py | 4 +- src/sparging/correlations.py | 14 +++---- src/sparging/inputs.py | 17 ++++----- src/sparging/model.py | 71 +++++++++++++++++++++++------------- test/test_solve.py | 53 ++++++++++++++++++++++++++- 6 files changed, 126 insertions(+), 52 deletions(-) diff --git a/sparging101.py b/sparging101.py index 2f13823..3608656 100644 --- a/sparging101.py +++ b/sparging101.py @@ -22,7 +22,7 @@ geom = ColumnGeometry( area=0.2 * ureg.m**2, - height=2 * ureg.m, + height=1 * ureg.m, nozzle_diameter=0.001 * ureg.m, nb_nozzle=10 * ureg.dimensionless, ) @@ -49,8 +49,9 @@ ) logger.info(my_input) -print(my_input.source_T_norm) -print(my_input.source_T_int.to("molT/s")) +print(my_input.get_S_T()) +print(f"{my_input.Q_T.to('molT/s')} = {my_input.Q_T.to('molT2/hour')}") +print(my_input.volume.to("m^3")) print( f"Concentration at steady state (no dispersion, no PP limited): { my_input.get_c_T2_SS().to('molT2/m^3') @@ -65,8 +66,10 @@ def profile_source_T(z: pint.Quantity | list[float], height: pint.Quantity = Non import numpy as np if isinstance(z, (float, np.ndarray, list)): # non-dimensional height (0 to 1) - return np.pi / 2 * np.sin(np.pi * z) # normalized + # return np.pi / 2 * np.sin(np.pi * z) # normalized + return 1 + 1 * np.sin(np.pi * z) # not normalized if isinstance(z, ureg.Quantity): + assert False if height is None: raise ValueError("Must provide height if z is a dimensional quantity") return np.pi / 2 * np.sin(np.pi / height * z) # normalized @@ -75,19 +78,21 @@ def profile_source_T(z: pint.Quantity | list[float], height: pint.Quantity = Non # return 0.5 * (1 + np.cos(0.5 * np.pi / (1 * ureg.m) * z)) +T_99 = 1 * ureg.hour my_simulation = Simulation( my_input, t_final=2 * T_99, signal_irr=lambda t: 1 if t < T_99 else 0, - signal_sparging=lambda t: 1, + signal_sparging=lambda t: 0, profile_pressure_hydrostatic=False, - # profile_source_T=profile_source_T, + profile_source_T=profile_source_T, + dispersion_on=False, ) if __name__ == "__main__": # my_simulation.sim_input.E_g *= 1e5 # my_simulation.sim_input.E_l *= 1e-5 - output = my_simulation.solve() + output = my_simulation.solve(fast_solve=True) # # save output to file # output.profiles_to_csv(f"output_{tank_height}m.csv") diff --git a/sparging_LIBRA1L.py b/sparging_LIBRA1L.py index 0fa2c77..fe29bb7 100644 --- a/sparging_LIBRA1L.py +++ b/sparging_LIBRA1L.py @@ -46,8 +46,8 @@ geom, flibe, operating_params, sparging_params ) logger.info(my_input) -print(my_input.source_T_norm) -print(my_input.source_T_int.to("molT/s")) +print(my_input.get_S_T()) +print(my_input.Q_T.to("molT/s")) n_fluence = 2.5e13 * ureg("neutron") n_gen_rate = operating_params.n_gen_rate diff --git a/src/sparging/correlations.py b/src/sparging/correlations.py index 7d41f32..e8f4c0c 100644 --- a/src/sparging/correlations.py +++ b/src/sparging/correlations.py @@ -393,16 +393,16 @@ def __contains__(self, key: str | Correlation): ) all_correlations.append(specific_interfacial_area) -source_T_normalized = Correlation( # TODO remettre source_T - identifier="source_T_norm", - function=lambda tbr, n_gen_rate, tank_volume: ( - tbr * n_gen_rate / tank_volume +source_T_integral = Correlation( + identifier="Q_T", + function=lambda tbr, n_gen_rate: ( + tbr * n_gen_rate ), # source term for tritium generation calculated from TBR and neutron generation rate corr_type=CorrelationType.TRITIUM_SOURCE, - input_units=["triton/neutron", "neutron/s", "m**3"], - output_units="molT/m**3/s", + input_units=["triton/neutron", "neutron/s"], + output_units="molT/s", ) -all_correlations.append(source_T_normalized) +all_correlations.append(source_T_integral) def get_d_b( diff --git a/src/sparging/inputs.py b/src/sparging/inputs.py index 29efe89..2b0aa9a 100644 --- a/src/sparging/inputs.py +++ b/src/sparging/inputs.py @@ -47,7 +47,7 @@ class OperatingParameters: P_bottom: pint.Quantity | Correlation | None = None tbr: pint.Quantity | None = None n_gen_rate: pint.Quantity | None = None - source_T_norm: pint.Quantity | Correlation | None = None + Q_T: pint.Quantity | None = None @dataclass @@ -77,7 +77,7 @@ class SimulationInput: E_g: pint.Quantity E_l: pint.Quantity D_l: pint.Quantity - source_T_norm: pint.Quantity + Q_T: pint.Quantity # normalize_source_T: bool = True @property @@ -88,21 +88,18 @@ def volume(self): def eps_l(self): return 1 - self.eps_g - @property # TODO enlever, inutile - def source_T_int(self): # rename Q_T2 - return self.source_T_norm * self.volume + def set_S_T(self, val: pint.Quantity): + self.Q_T = (val.to("molT/m**3/s") * self.volume).to("molT/s") - def set_S_T2(self, val: pint.Quantity): # TODO to use - self.Q_T2 = (val.to("molT2/m**3/s") * self.volume).to("molT2/s") + def get_S_T(self) -> pint.Quantity: + return (self.Q_T / self.volume).to("molT/m**3/s") def get_tau(self) -> pint.Quantity: """characteristic time of the sparger under the small partial pressure (SPP) approximation""" return (self.eps_l / (self.h_l * self.a)).to("seconds") def get_c_T2_SS(self) -> pint.Quantity: - return (self.source_T_int / self.volume * 1 / (self.h_l * self.a)).to( - "molT2/m^3" - ) + return (self.get_S_T() * 1 / (self.h_l * self.a)).to("molT2/m^3") def get_PP_number(self) -> pint.Quantity: """Partial pressure number, ratio of the equivalent T concentration at liquid boundary to the bulk liquid concentration. diff --git a/src/sparging/model.py b/src/sparging/model.py index 8ae15cc..ed94606 100644 --- a/src/sparging/model.py +++ b/src/sparging/model.py @@ -143,7 +143,9 @@ class Simulation: signal_irr: callable[pint.Quantity] signal_sparging: callable[pint.Quantity] profile_source_T: callable[pint.Quantity] | None = None + """callable = f:[0,1] -> R+, it takes a dimensionless coordinate: (z / height)""" profile_pressure_hydrostatic: bool = False + dispersion_on: bool = True def hydrostatic_pressure(self, x: pint.Quantity) -> pint.Quantity: """returns the hydrostatic pressure at a given height x in the tank given P_bottom""" @@ -151,7 +153,12 @@ def hydrostatic_pressure(self, x: pint.Quantity) -> pint.Quantity: g = const_g return (self.sim_input.P_bottom + rho * g * x).to("Pa") - def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None): + def solve( + self, + dt: pint.Quantity | None = None, + dx: pint.Quantity | None = None, + fast_solve: bool = False, + ) -> SimulationResults: # unpack pint.Quantities t_final = self.t_final.to("seconds").magnitude tank_height = self.sim_input.height.to("m").magnitude @@ -167,10 +174,18 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None E_l = self.sim_input.E_l.to("m**2/s").magnitude D_l = self.sim_input.D_l.to("m**2/s").magnitude # not needed (included in h_l) u_g0 = self.sim_input.u_g0.to("m/s").magnitude - source_T2_norm = self.sim_input.source_T_norm.to("molT2/s/m**3").magnitude + Q_T2 = self.sim_input.Q_T.to("molT2/s").magnitude - dt = dt.to("seconds").magnitude if dt is not None else t_final / 1000 - dx = dx.to("m").magnitude if dx is not None else tank_height / 1000 + dt = ( + dt.to("seconds").magnitude + if dt is not None + else (t_final / 1000 if not fast_solve else t_final / 50) + ) + dx = ( + dx.to("m").magnitude + if dx is not None + else (tank_height / 1000 if not fast_solve else tank_height / 50) + ) eps_l = 1 - eps_g # MESH AND FUNCTION SPACES @@ -196,23 +211,27 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None h_l_const = dolfinx.fem.Constant(mesh, PETSc.ScalarType(h_l)) - gen_T2_norm = dolfinx.fem.Constant( - mesh, source_T2_norm * self.signal_irr(0 * ureg.s) + gen_T2_ave = dolfinx.fem.Constant( + mesh, Q_T2 / tank_volume * self.signal_irr(0 * ureg.s) ) # magnitude of the generation term if self.profile_source_T is not None: # spatially varying profile is provided - gen_T2_prof = dolfinx.fem.Function(V_profile) - gen_T2_prof.interpolate( + arbitrary_profile = dolfinx.fem.Function(V_profile) + arbitrary_profile.interpolate( lambda x: x[0] * 0 + self.profile_source_T(x[0] / tank_height) ) + profile_integral = dolfinx.fem.assemble_scalar( + dolfinx.fem.form( + arbitrary_profile * 1 / tank_height * ufl.dx # TODO + ) # dimensionless integral + ) + normalized_profile = ( + arbitrary_profile / profile_integral + ) # normalize profile so that its integral over the dimensionless height is 1 else: # homogeneous generation - gen_T2_prof = dolfinx.fem.Constant(mesh, PETSc.ScalarType(1.0)) + normalized_profile = dolfinx.fem.Constant(mesh, PETSc.ScalarType(1.0)) - # if self.sim_input.normalize_source_T: - gen_T2 = ( - gen_T2_norm * gen_T2_prof - ) # gen_T2_norm is already mormalized by tank volume - # user can choose to specify a normalized profile_source_T so that mean(gen_T2) = gen_T2_norm = 0.5 * source_T_norm + gen_T2 = gen_T2_ave * normalized_profile P_prof = dolfinx.fem.Function(V_profile) if self.profile_pressure_hydrostatic: @@ -238,15 +257,16 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None F += eps_g * 1 / (const.R * T) * (P * (y_T2 - y_T2_n) / dt) * v_y * ufl.dx # dispersive terms - F += eps_l * E_l * ufl.dot(ufl.grad(c_T2), ufl.grad(v_c)) * ufl.dx - F += ( - eps_g - * E_g - * 1 - / (const.R * T) - * ufl.dot(ufl.grad(P * y_T2), ufl.grad(v_y)) - * ufl.dx - ) + if self.dispersion_on is True: + F += eps_l * E_l * ufl.dot(ufl.grad(c_T2), ufl.grad(v_c)) * ufl.dx + F += ( + eps_g + * E_g + * 1 + / (const.R * T) + * ufl.dot(ufl.grad(P * y_T2), ufl.grad(v_y)) + * ufl.dx + ) # mass exchange (coupling term) F += J_T2 * v_c * ufl.dx - J_T2 * v_y * ufl.dx @@ -309,6 +329,7 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None y_sort_coords = np.argsort(coords) x_y = coords[y_sort_coords] + # TODO Initialize results storage with zeros -> there's an inconsistency: times[0] = 0 but inventories[0] != 0 -> screws comparison with analytical solutions times = [] c_T2_solutions = [] y_T2_solutions = [] @@ -321,7 +342,7 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None t = 0 while t < t_final: # update time-dependent terms - gen_T2_norm.value = source_T2_norm * self.signal_irr(t * ureg.s) + gen_T2_ave.value = Q_T2 / tank_volume * self.signal_irr(t * ureg.s) h_l_const.value = h_l * self.signal_sparging(t * ureg.s) problem.solve() @@ -345,7 +366,7 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None y_T2_solutions.append(y_T2_vals.copy()) J_T2_solutions.append(J_T2_vals.copy()) sources_T2.append( - source_T2_norm * self.signal_irr(t * ureg.s) * tank_volume + Q_T2 * self.signal_irr(t * ureg.s) ) # total T generation rate in the tank [mol/s] TODO useless: signal_irr is already given n = ufl.FacetNormal(mesh) diff --git a/test/test_solve.py b/test/test_solve.py index 2894840..963fdd4 100644 --- a/test/test_solve.py +++ b/test/test_solve.py @@ -4,6 +4,7 @@ import pytest import dataclasses from pint import DimensionalityError +import numpy as np standard_input = SimulationInput( @@ -20,7 +21,7 @@ E_g=1e-2 * ureg("m^2/s"), E_l=1e-1 * ureg("m^2/s"), D_l=3e-9 * ureg("m^2/s"), - source_T=8e-16 * ureg("molT/m^3/s"), + Q_T=1e8 * ureg("T/s"), ) @@ -85,3 +86,53 @@ def test_model_solve_wrong_argument(standard_simulation): """ with pytest.raises(AttributeError, match="object has no attribute 'to'"): standard_simulation.solve(dt=0.01) + + +@pytest.mark.parametrize("case", ("nominal", "volume_change", "unormalized_profile")) +def test_source_T_normalization(case): + """Tests that the source term is correctly normalized so that the total inventory after a fixed time is always the same + no matter the volume or the given profile function.""" + Q_T = standard_input.Q_T + t_final = 50 * ureg.seconds + # BUILD + match case: + case "nominal": + my_simulation = Simulation( + standard_input, + t_final=t_final, + signal_irr=lambda t: 1, + signal_sparging=lambda t: 0, + ) + case "volume_change": + modified_input = dataclasses.replace( + standard_input, + height=23.2 * standard_input.height, + area=0.123 * standard_input.area, + ) + my_simulation = Simulation( + modified_input, + t_final=t_final, + signal_irr=lambda t: 1, + signal_sparging=lambda t: 0, + ) + case "unormalized_profile": + my_simulation = Simulation( + standard_input, + t_final=t_final, + signal_irr=lambda t: 1, + signal_sparging=lambda t: 0, + profile_source_T=lambda xi: 3 + xi, # not normalized + ) + # RUN + output = my_simulation.solve(fast_solve=True) + + # TEST + n_result = ureg.Quantity(output.inventories_T2_salt[-1], "molT2").magnitude + # n_theory = (Q_T * ureg.Quantity(output.times[-1], "s")).to("molT2").magnitude # TODO change to this + n_theory = (Q_T * t_final).to("molT2").magnitude + # print(len(output.times), len(output.inventories_T2_salt)) + # print(output.times[0], output.inventories_T2_salt[0]) + # print(output.times[-1], output.inventories_T2_salt[-1]) + assert np.isclose(n_result, n_theory, atol=0, rtol=1e-2), print( + f"n_result = {n_result}, n_theory = {n_theory}" + ) # TODO reduce tolerance when vectors initialization with zero is fixed From c3ea7413bb642566c133623b75511812cd56b267 Mon Sep 17 00:00:00 2001 From: hmonroedd Date: Thu, 16 Apr 2026 12:09:08 -0400 Subject: [PATCH 10/12] refactoring with postprocess function --- src/sparging/model.py | 56 +++++++++++++++++++------------------------ 1 file changed, 25 insertions(+), 31 deletions(-) diff --git a/src/sparging/model.py b/src/sparging/model.py index ed94606..1e2a0d9 100644 --- a/src/sparging/model.py +++ b/src/sparging/model.py @@ -294,11 +294,6 @@ def solve( dolfinx.fem.locate_dofs_topological(V.sub(1), fdim, gas_inlet_facets), V.sub(1), ) # y_T2 = 0 at gas inlet - # bc2 = dolfinx.fem.dirichletbc( - # dolfinx.fem.Constant(mesh, 0.0), - # dolfinx.fem.locate_dofs_topological(V.sub(0), fdim, gas_outlet_facets), - # V.sub(0), - # ) # c_T2 = 0 at gas outlet # Custom measure all_facets = np.concatenate((gas_inlet_facets, gas_outlet_facets)) @@ -312,7 +307,6 @@ def solve( problem = NonlinearProblem( F, u, - # bcs=[bc1, bc2], bcs=[bc1], # Neumann BCs on c_T2 at inlet and outlet are naturally enforced petsc_options_prefix="librasparge", # petsc_options={"snes_monitor": None}, @@ -329,28 +323,7 @@ def solve( y_sort_coords = np.argsort(coords) x_y = coords[y_sort_coords] - # TODO Initialize results storage with zeros -> there's an inconsistency: times[0] = 0 but inventories[0] != 0 -> screws comparison with analytical solutions - times = [] - c_T2_solutions = [] - y_T2_solutions = [] - J_T2_solutions = [] - sources_T2 = [] - fluxes_T2 = [] - inventories_T2_salt = [] - - # SOLVE - t = 0 - while t < t_final: - # update time-dependent terms - gen_T2_ave.value = Q_T2 / tank_volume * self.signal_irr(t * ureg.s) - h_l_const.value = h_l * self.signal_sparging(t * ureg.s) - - problem.solve() - - # update previous solution - u_n.x.array[:] = u.x.array[:] - - # post process + def post_process(t): c_T2_post, y_T2_post = u.split() c_T2_vals = u.x.array[ct_dofs][ct_sort_coords] @@ -358,9 +331,6 @@ def solve( J_T2_vals = ( a * h_l_const.value * (c_T2_vals - K_s * (P.x.array * y_T2_vals + EPS)) ) # TODO there is a clever way of getting J_T2 for sure - - # store time and solution - # TODO give units to the results times.append(t) c_T2_solutions.append(c_T2_vals.copy()) y_T2_solutions.append(y_T2_vals.copy()) @@ -410,6 +380,30 @@ def solve( inventory_T2_salt *= tank_area # get total amount of T2 in [mol] inventories_T2_salt.append(inventory_T2_salt) + # TODO Initialize results storage with zeros -> there's an inconsistency: times[0] = 0 but inventories[0] != 0 -> screws comparison with analytical solutions + t = 0 + times = [] + c_T2_solutions = [] + y_T2_solutions = [] + J_T2_solutions = [] + sources_T2 = [] + fluxes_T2 = [] + inventories_T2_salt = [] + post_process(t) + + # SOLVE + while t < t_final: + # update time-dependent terms + gen_T2_ave.value = Q_T2 / tank_volume * self.signal_irr(t * ureg.s) + h_l_const.value = h_l * self.signal_sparging(t * ureg.s) + + problem.solve() + + # update previous solution + u_n.x.array[:] = u.x.array[:] + + post_process(t) + # advance time t += dt From 5edb024e5e6d2710281cf862d52ab2e072ffa7a7 Mon Sep 17 00:00:00 2001 From: hmonroedd Date: Thu, 16 Apr 2026 18:08:49 -0400 Subject: [PATCH 11/12] cleaner J_T2 post processing + retrocompatibility with python 3.12 to use with autoemulate --- autoemulate_env.yml | 15 +++++++++++++++ sparging101.py | 8 ++++---- src/sparging/config.py | 2 +- src/sparging/inputs.py | 1 + src/sparging/model.py | 33 +++++++++++++++++++++++++++------ 5 files changed, 48 insertions(+), 11 deletions(-) create mode 100644 autoemulate_env.yml diff --git a/autoemulate_env.yml b/autoemulate_env.yml new file mode 100644 index 0000000..f2c07c9 --- /dev/null +++ b/autoemulate_env.yml @@ -0,0 +1,15 @@ +name: autoemulate_env +channels: + - conda-forge +dependencies: + - python=3.12 + - fenics-dolfinx + - matplotlib + - pyvista + - pyyaml + - numpy + - scipy + - pandas + - pint + - pip: + - autoemulate \ No newline at end of file diff --git a/sparging101.py b/sparging101.py index 3608656..231559e 100644 --- a/sparging101.py +++ b/sparging101.py @@ -1,3 +1,4 @@ +from __future__ import annotations from sparging.config import ureg from sparging import all_correlations from sparging import animation @@ -78,21 +79,20 @@ def profile_source_T(z: pint.Quantity | list[float], height: pint.Quantity = Non # return 0.5 * (1 + np.cos(0.5 * np.pi / (1 * ureg.m) * z)) -T_99 = 1 * ureg.hour +T_99 = my_input.get_tau() my_simulation = Simulation( my_input, t_final=2 * T_99, signal_irr=lambda t: 1 if t < T_99 else 0, - signal_sparging=lambda t: 0, + signal_sparging=lambda t: 1, profile_pressure_hydrostatic=False, profile_source_T=profile_source_T, - dispersion_on=False, ) if __name__ == "__main__": # my_simulation.sim_input.E_g *= 1e5 # my_simulation.sim_input.E_l *= 1e-5 - output = my_simulation.solve(fast_solve=True) + output = my_simulation.solve() # # save output to file # output.profiles_to_csv(f"output_{tank_height}m.csv") diff --git a/src/sparging/config.py b/src/sparging/config.py index 1ec8f5c..bf9bd10 100644 --- a/src/sparging/config.py +++ b/src/sparging/config.py @@ -9,7 +9,7 @@ ureg.define(f"molT = {const.N_A} * triton") ureg.define(f"molT2 = 2 * {const.N_A} * triton") ureg.define("neutron = [neutron] = n") -ureg.define("sccm = 7.44e-7 mol/s") +ureg.define("sccm = 7.44e-7 mol/s") # holds for an ideal gas const_R = const.R * ureg("J/K/mol") # ideal gas constant diff --git a/src/sparging/inputs.py b/src/sparging/inputs.py index 2b0aa9a..28a5597 100644 --- a/src/sparging/inputs.py +++ b/src/sparging/inputs.py @@ -79,6 +79,7 @@ class SimulationInput: D_l: pint.Quantity Q_T: pint.Quantity # normalize_source_T: bool = True + # TODO refactoring, put signals in this class @property def volume(self): diff --git a/src/sparging/model.py b/src/sparging/model.py index 1e2a0d9..49ca119 100644 --- a/src/sparging/model.py +++ b/src/sparging/model.py @@ -1,3 +1,4 @@ +from __future__ import annotations from mpi4py import MPI import dolfinx import basix @@ -24,6 +25,7 @@ if TYPE_CHECKING: import pint +from collections.abc import Callable hours_to_seconds = 3600 days_to_seconds = 24 * hours_to_seconds @@ -140,9 +142,9 @@ def profiles_to_csv(self, output_path: Path): class Simulation: sim_input: SimulationInput t_final: pint.Quantity - signal_irr: callable[pint.Quantity] - signal_sparging: callable[pint.Quantity] - profile_source_T: callable[pint.Quantity] | None = None + signal_irr: Callable[[pint.Quantity], float] + signal_sparging: Callable[[pint.Quantity], float] + profile_source_T: Callable[[float], float] | None = None """callable = f:[0,1] -> R+, it takes a dimensionless coordinate: (z / height)""" profile_pressure_hydrostatic: bool = False dispersion_on: bool = True @@ -313,6 +315,11 @@ def solve( ) # initialise post processing + + # we define a function and an expression for J_T2 for use in post processing + J_T2_func = dolfinx.fem.Function(V_profile) + J_T2_expr = dolfinx.fem.Expression(J_T2, V_profile.element.interpolation_points) + V0_ct, ct_dofs = u.function_space.sub(0).collapse() coords = V0_ct.tabulate_dof_coordinates()[:, 0] ct_sort_coords = np.argsort(coords) @@ -323,14 +330,28 @@ def solve( y_sort_coords = np.argsort(coords) x_y = coords[y_sort_coords] + coords_profile = V_profile.tabulate_dof_coordinates()[:, 0] + num_profile_dofs = ( + V_profile.dofmap.index_map.size_local + ) * V_profile.dofmap.index_map_bs + + profile_dofs = np.arange(num_profile_dofs) + profile_sort_coords = np.argsort(coords_profile) + # NOTE currently we don't use x_profile and use another x in the plotting script + x_profile = coords_profile[profile_sort_coords] + + # NOTE maybe we could take this function out and it would take a SimulationResults object as input + u + other things... def post_process(t): + """ + Post-process the solution at time t. + Extract solution profiles, compute fluxes and inventories, and store them in lists for later analysis. + """ c_T2_post, y_T2_post = u.split() c_T2_vals = u.x.array[ct_dofs][ct_sort_coords] y_T2_vals = u.x.array[y_dofs][y_sort_coords] - J_T2_vals = ( - a * h_l_const.value * (c_T2_vals - K_s * (P.x.array * y_T2_vals + EPS)) - ) # TODO there is a clever way of getting J_T2 for sure + J_T2_func.interpolate(J_T2_expr) + J_T2_vals = J_T2_func.x.array[profile_dofs][ct_sort_coords] times.append(t) c_T2_solutions.append(c_T2_vals.copy()) y_T2_solutions.append(y_T2_vals.copy()) From b66ec0518d04e3ad85be76efe07d3a66b2f7cbfb Mon Sep 17 00:00:00 2001 From: hmonroedd Date: Fri, 17 Apr 2026 16:31:59 -0400 Subject: [PATCH 12/12] removed unused file --- autoemulate_env.yml | 15 --------------- 1 file changed, 15 deletions(-) delete mode 100644 autoemulate_env.yml diff --git a/autoemulate_env.yml b/autoemulate_env.yml deleted file mode 100644 index f2c07c9..0000000 --- a/autoemulate_env.yml +++ /dev/null @@ -1,15 +0,0 @@ -name: autoemulate_env -channels: - - conda-forge -dependencies: - - python=3.12 - - fenics-dolfinx - - matplotlib - - pyvista - - pyyaml - - numpy - - scipy - - pandas - - pint - - pip: - - autoemulate \ No newline at end of file