From 04ad1aa077be6a28565eec82f686f0db00983143 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 19 Dec 2024 17:43:39 -0500 Subject: [PATCH 01/11] added a test to catch the bug --- test/test_helpers.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) create mode 100644 test/test_helpers.py diff --git a/test/test_helpers.py b/test/test_helpers.py new file mode 100644 index 0000000..d71410d --- /dev/null +++ b/test/test_helpers.py @@ -0,0 +1,12 @@ +from hisp.helpers import gaussian_distribution + +import numpy as np + + +def test_gaussian_distr_area(): + x = np.linspace(-10, 10, 1000) + y = gaussian_distribution([x], mean=0, width=1, mod=np) + area = np.trapz(y, x) + assert np.isclose( + area, 1, atol=1e-2 + ), f"Area under Gaussian distribution is {area} instead of 1" From 7e9f73a797b2f6f8404ea8e394b9e3de85f375c6 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 19 Dec 2024 17:46:54 -0500 Subject: [PATCH 02/11] fixed bug --- src/hisp/helpers.py | 4 +++- test/test_helpers.py | 1 + 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/hisp/helpers.py b/src/hisp/helpers.py index e6ad6fd..8a3acab 100644 --- a/src/hisp/helpers.py +++ b/src/hisp/helpers.py @@ -40,7 +40,9 @@ def update(self, t: float): def gaussian_distribution( x: npt.NDArray, mean: float, width: float, mod=ufl ) -> ufl.core.expr.Expr: - return mod.exp(-((x[0] - mean) ** 2) / (2 * width**2)) + return mod.exp(-((x[0] - mean) ** 2) / (2 * width**2)) / ( + np.sqrt(2 * np.pi * width**2) + ) def periodic_step_function(x, period_on, period_total, value, value_off=0.0): diff --git a/test/test_helpers.py b/test/test_helpers.py index d71410d..92eb91b 100644 --- a/test/test_helpers.py +++ b/test/test_helpers.py @@ -4,6 +4,7 @@ def test_gaussian_distr_area(): + """Generates a gaussian distribution and checks if the area under the curve is 1.""" x = np.linspace(-10, 10, 1000) y = gaussian_distribution([x], mean=0, width=1, mod=np) area = np.trapz(y, x) From 52b7e9f44580aa7066aa5f605aa5acf0600235ee Mon Sep 17 00:00:00 2001 From: kaelyndunnell Date: Thu, 19 Dec 2024 21:18:33 -0800 Subject: [PATCH 03/11] add photon radiation heat to plasmadatahandling class (initialized at a given value) --- src/hisp/plamsa_data_handling/main.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/hisp/plamsa_data_handling/main.py b/src/hisp/plamsa_data_handling/main.py index 924febe..a2aeb77 100644 --- a/src/hisp/plamsa_data_handling/main.py +++ b/src/hisp/plamsa_data_handling/main.py @@ -223,7 +223,7 @@ def get_heat(self, pulse: Pulse, bin: SubBin | DivBin, t_rel: float) -> float: raise ValueError(f"Invalid pulse type {pulse.pulse_type}") if pulse.pulse_type == "FP": - heat_total = data["heat_total"][bin_index] + heat_total = data["heat_total"][bin_index]+0.11e+06 # photon radiation heat adjustment heat_ion = data["heat_ion"][bin_index] if isinstance(bin, SubBin): heat_val = heat_total - heat_ion * (1 - bin.wetted_frac) @@ -231,7 +231,7 @@ def get_heat(self, pulse: Pulse, bin: SubBin | DivBin, t_rel: float) -> float: heat_val = heat_total elif pulse.pulse_type == "RISP": if isinstance(bin, SubBin): - heat_total = data["heat_total"] + heat_total = data["heat_total"]+0.11e+06 # photon radiation heat adjustment heat_ion = data["heat_ion"] heat_val = heat_total - heat_ion * (1 - bin.wetted_frac) else: From 446a07891e87e09b89c450d861a8b4697ef36e2e Mon Sep 17 00:00:00 2001 From: kaelyndunnell Date: Thu, 19 Dec 2024 21:25:20 -0800 Subject: [PATCH 04/11] realistic tolerances --- src/hisp/festim_models/mb_model.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/hisp/festim_models/mb_model.py b/src/hisp/festim_models/mb_model.py index d79eeab..cc00f5b 100644 --- a/src/hisp/festim_models/mb_model.py +++ b/src/hisp/festim_models/mb_model.py @@ -274,7 +274,7 @@ def make_W_mb_model( ############# Settings ############# my_model.settings = F.Settings( - atol=1e-10, + atol=1e10, rtol=1e-10, max_iterations=1000, final_time=final_time, @@ -538,8 +538,8 @@ def make_B_mb_model( ############# Settings ############# my_model.settings = F.Settings( - atol=1e-15, - rtol=1e-15, + atol=1e10, + rtol=1e-10, max_iterations=1000, final_time=final_time, ) @@ -743,8 +743,8 @@ def make_DFW_mb_model( ############# Settings ############# my_model.settings = F.Settings( - atol=1e-15, - rtol=1e-15, + atol=1e10, + rtol=1e-10, max_iterations=1000, final_time=final_time, ) From f46e7f9b5a4e7845e5bc2277f8e25a7340bb1dac Mon Sep 17 00:00:00 2001 From: kaelyndunnell Date: Thu, 19 Dec 2024 21:40:31 -0800 Subject: [PATCH 05/11] W mesh refinement --- src/hisp/festim_models/mb_model.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/hisp/festim_models/mb_model.py b/src/hisp/festim_models/mb_model.py index cc00f5b..5887b1d 100644 --- a/src/hisp/festim_models/mb_model.py +++ b/src/hisp/festim_models/mb_model.py @@ -48,10 +48,11 @@ def make_W_mb_model( vertices = np.concatenate( # 1D mesh with extra refinement [ - np.linspace(0, 30e-9, num=200), - np.linspace(30e-9, 3e-6, num=300), - np.linspace(3e-6, 30e-6, num=200), - np.linspace(30e-6, L, num=200), + np.linspace(0, 30e-9, num=300), + np.linspace(30e-9, 3e-6, num=400), + np.linspace(3e-6, 30e-6, num=400), + np.linspace(30e-6, 1e-4, num=400), + np.linspace(1e-4, L, num=300), ] ) my_model.mesh = F.Mesh1D(vertices) From 8ffb867bb4114ab22f4bad6bb8cb3875f37b31e0 Mon Sep 17 00:00:00 2001 From: kaelyndunnell Date: Thu, 19 Dec 2024 22:08:24 -0800 Subject: [PATCH 06/11] updated boron mesh and tolerances --- src/hisp/festim_models/mb_model.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/hisp/festim_models/mb_model.py b/src/hisp/festim_models/mb_model.py index 5887b1d..800bfa4 100644 --- a/src/hisp/festim_models/mb_model.py +++ b/src/hisp/festim_models/mb_model.py @@ -318,8 +318,9 @@ def make_B_mb_model( vertices = np.concatenate( # 1D mesh with extra refinement [ - np.linspace(0, 30e-9, num=200), - np.linspace(30e-9, L, num=200), + np.linspace(0, 5e-9, num=300), + np.linspace(5e-9, 1e-8, num=500), + np.linspace(1e-8, L, num=300), ] ) my_model.mesh = F.Mesh1D(vertices) @@ -540,7 +541,7 @@ def make_B_mb_model( ############# Settings ############# my_model.settings = F.Settings( atol=1e10, - rtol=1e-10, + rtol=1e-5, max_iterations=1000, final_time=final_time, ) From 3d6655a704ad8f1ec1cc51248a617aefbb248fa1 Mon Sep 17 00:00:00 2001 From: Kaelyn Dunnell <150201080+kaelyndunnell@users.noreply.github.com> Date: Fri, 20 Dec 2024 15:11:42 -0800 Subject: [PATCH 07/11] Update src/hisp/plamsa_data_handling/main.py MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Rémi Delaporte-Mathurin <40028739+RemDelaporteMathurin@users.noreply.github.com> --- src/hisp/plamsa_data_handling/main.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/hisp/plamsa_data_handling/main.py b/src/hisp/plamsa_data_handling/main.py index a2aeb77..6621227 100644 --- a/src/hisp/plamsa_data_handling/main.py +++ b/src/hisp/plamsa_data_handling/main.py @@ -223,7 +223,8 @@ def get_heat(self, pulse: Pulse, bin: SubBin | DivBin, t_rel: float) -> float: raise ValueError(f"Invalid pulse type {pulse.pulse_type}") if pulse.pulse_type == "FP": - heat_total = data["heat_total"][bin_index]+0.11e+06 # photon radiation heat adjustment + photon_heat_radiation = 0.11e6 # W/m2 + heat_total = data["heat_total"][bin_index]+ photon_heat_radiation heat_ion = data["heat_ion"][bin_index] if isinstance(bin, SubBin): heat_val = heat_total - heat_ion * (1 - bin.wetted_frac) From b247f434b8e9c1b765f668d33f3a824ae06dbb56 Mon Sep 17 00:00:00 2001 From: Kaelyn Dunnell <150201080+kaelyndunnell@users.noreply.github.com> Date: Fri, 20 Dec 2024 15:12:12 -0800 Subject: [PATCH 08/11] Update src/hisp/plamsa_data_handling/main.py MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Rémi Delaporte-Mathurin <40028739+RemDelaporteMathurin@users.noreply.github.com> --- src/hisp/plamsa_data_handling/main.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/hisp/plamsa_data_handling/main.py b/src/hisp/plamsa_data_handling/main.py index 6621227..e766e4e 100644 --- a/src/hisp/plamsa_data_handling/main.py +++ b/src/hisp/plamsa_data_handling/main.py @@ -232,7 +232,8 @@ def get_heat(self, pulse: Pulse, bin: SubBin | DivBin, t_rel: float) -> float: heat_val = heat_total elif pulse.pulse_type == "RISP": if isinstance(bin, SubBin): - heat_total = data["heat_total"]+0.11e+06 # photon radiation heat adjustment + photon_radiation_heat = 0.11e6 # W/m2 + heat_total = data["heat_total"]+photon_radiation_heat heat_ion = data["heat_ion"] heat_val = heat_total - heat_ion * (1 - bin.wetted_frac) else: From e0b1eed661e8494b97be2bc47797cebcc17fc416 Mon Sep 17 00:00:00 2001 From: kaelyndunnell Date: Fri, 20 Dec 2024 15:19:23 -0800 Subject: [PATCH 09/11] added documentation to gaussian_distribution --- src/hisp/helpers.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/hisp/helpers.py b/src/hisp/helpers.py index 8a3acab..894109c 100644 --- a/src/hisp/helpers.py +++ b/src/hisp/helpers.py @@ -40,6 +40,17 @@ def update(self, t: float): def gaussian_distribution( x: npt.NDArray, mean: float, width: float, mod=ufl ) -> ufl.core.expr.Expr: + """Generates a gaussian distribution for particle sources. + + Args: + x (npt.NDArray): x values along the length of given bin. + mean (float): Mean of the distribution. + width (float): Width of the gaussian distribution. + mod (_type_, optional): Module used to express gaussian distribution. Defaults to ufl. + + Returns: + ufl.core.expr.Expr: Gaussian distribution with area 1. + """ return mod.exp(-((x[0] - mean) ** 2) / (2 * width**2)) / ( np.sqrt(2 * np.pi * width**2) ) From 04a2e2c9b344699f391c84c14370c9a7771557c0 Mon Sep 17 00:00:00 2001 From: kaelyndunnell Date: Fri, 20 Dec 2024 15:21:30 -0800 Subject: [PATCH 10/11] decrease some B mesh refinement --- src/hisp/festim_models/mb_model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/hisp/festim_models/mb_model.py b/src/hisp/festim_models/mb_model.py index 800bfa4..b77d724 100644 --- a/src/hisp/festim_models/mb_model.py +++ b/src/hisp/festim_models/mb_model.py @@ -318,7 +318,7 @@ def make_B_mb_model( vertices = np.concatenate( # 1D mesh with extra refinement [ - np.linspace(0, 5e-9, num=300), + np.linspace(0, 5e-9, num=200), np.linspace(5e-9, 1e-8, num=500), np.linspace(1e-8, L, num=300), ] From 1d282b87f97ad3e5dbe0ebc274f7d802270c6875 Mon Sep 17 00:00:00 2001 From: kaelyndunnell Date: Fri, 20 Dec 2024 15:25:34 -0800 Subject: [PATCH 11/11] decrease B model rtol --- src/hisp/festim_models/mb_model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/hisp/festim_models/mb_model.py b/src/hisp/festim_models/mb_model.py index b77d724..870fe58 100644 --- a/src/hisp/festim_models/mb_model.py +++ b/src/hisp/festim_models/mb_model.py @@ -541,7 +541,7 @@ def make_B_mb_model( ############# Settings ############# my_model.settings = F.Settings( atol=1e10, - rtol=1e-5, + rtol=1e-8, max_iterations=1000, final_time=final_time, )