From 0e3004d0511aec56faa6356ba57e8b2598a6a2c9 Mon Sep 17 00:00:00 2001 From: KulaginVladimir Date: Mon, 9 Sep 2024 21:26:52 +0300 Subject: [PATCH 1/8] added a write_at_last parameter --- festim/__init__.py | 2 +- festim/exports/exports.py | 3 +- festim/exports/txt_export.py | 107 +++++++++------------ test/unit/test_exports/test_txt_export.py | 31 +----- test/unit/test_exports/test_txt_exports.py | 28 ------ 5 files changed, 52 insertions(+), 119 deletions(-) delete mode 100644 test/unit/test_exports/test_txt_exports.py diff --git a/festim/__init__.py b/festim/__init__.py index 10ce0105c..1609104cb 100644 --- a/festim/__init__.py +++ b/festim/__init__.py @@ -88,7 +88,7 @@ from .exports.derived_quantities.derived_quantities import DerivedQuantities -from .exports.txt_export import TXTExport, TXTExports +from .exports.txt_export import TXTExport from .settings import Settings diff --git a/festim/exports/exports.py b/festim/exports/exports.py index 45e769a62..c3d61480a 100644 --- a/festim/exports/exports.py +++ b/festim/exports/exports.py @@ -129,8 +129,7 @@ def write(self, label_to_function, dx): label_to_function[export.field], self.V_DG1 ) export.function = label_to_function[export.field] - steady = self.final_time == None - export.write(self.t, steady) + export.write(self.t, self.final_time) self.nb_iterations += 1 def initialise_derived_quantities(self, dx, ds, materials): diff --git a/festim/exports/txt_export.py b/festim/exports/txt_export.py index 3c7688a99..b306b31e3 100644 --- a/festim/exports/txt_export.py +++ b/festim/exports/txt_export.py @@ -20,16 +20,22 @@ class TXTExport(festim.Export): Defautls to ".2e". """ - def __init__(self, field, filename, times=None, header_format=".2e") -> None: + def __init__( + self, field, filename, times=None, write_at_last=False, header_format=".2e" + ) -> None: super().__init__(field=field) if times: self.times = sorted(times) else: self.times = times self.filename = filename + self.write_at_last = write_at_last self.header_format = header_format self._first_time = True + self.data = None + self.header = None + @property def filename(self): return self._filename @@ -51,15 +57,21 @@ def is_it_time_to_export(self, current_time): return True return False - def when_is_next_time(self, current_time): - if self.times is None: - return None - for time in self.times: - if current_time < time: - return time - return None + def is_last(self, current_time, final_time): + if final_time is None: + # write if steady + return True + elif self.times is None: + if np.isclose(current_time, final_time, atol=0): + # write at final time if exports at each timestep + return True + else: + if np.isclose(current_time, self.times[-1], atol=0): + # write at final time if exports at specific times + return True + return False - def write(self, current_time, steady): + def write(self, current_time, final_time): # create a DG1 functionspace V_DG1 = f.FunctionSpace(self.function.function_space().mesh(), "DG", 1) @@ -72,62 +84,35 @@ def write(self, current_time, steady): if not os.path.exists(dirname): os.makedirs(dirname, exist_ok=True) - # if steady or it is the first time to export - # write data + # write data if steady or it is the first time to export # else append new column to the existing file - if steady or self._first_time: - if steady: - header = "x,t=steady" + if final_time is None or self._first_time: + if final_time is None: + self.header = "x,t=steady" else: - header = f"x,t={format(current_time, self.header_format)}s" + self.header = f"x,t={format(current_time, self.header_format)}s" + x = f.interpolate(f.Expression("x[0]", degree=1), V_DG1) x_column = np.transpose([x.vector()[:]]) - data = np.column_stack([x_column, solution_column]) + self.data = np.column_stack([x_column, solution_column]) self._first_time = False else: # Update the header - old_file = open(self.filename) - old_header = old_file.readline().split("\n")[0] - old_file.close() - header = old_header + f",t={format(current_time, self.header_format)}s" - # Append new column - old_columns = np.loadtxt(self.filename, delimiter=",", skiprows=1) - data = np.column_stack([old_columns, solution_column]) - - np.savetxt(self.filename, data, header=header, delimiter=",", comments="") - - -class TXTExports: - """ - Args: - fields (list): list of exported fields ("solute", "1", "retention", - "T"...) - filenames (list): list of the filenames for each field (must end with .txt). - times (list, optional): if provided, fields will be - exported at these timesteps. Otherwise exports at all - timesteps. Defaults to None. - header_format (str, optional): the format of column headers. - Defautls to ".2e". - """ - - def __init__( - self, fields=[], filenames=[], times=None, header_format=".2e" - ) -> None: - msg = "TXTExports class will be deprecated in future versions of FESTIM" - warnings.warn(msg, DeprecationWarning) - - self.fields = fields - if len(self.fields) != len(filenames): - raise ValueError( - "Number of fields to be exported " - "doesn't match number of filenames in txt exports" - ) - if times: - self.times = sorted(times) - else: - self.times = times - self.filenames = filenames - self.header_format = header_format - self.exports = [] - for function, filename in zip(self.fields, self.filenames): - self.exports.append(TXTExport(function, filename, times, header_format)) + self.header += f",t={format(current_time, self.header_format)}s" + # Add new column + self.data = np.column_stack([self.data, solution_column]) + + if ( + self.write_at_last and self.is_last(current_time, final_time) + ) or not self.write_at_last: + if self.is_last(current_time, final_time): + # Sort data by the x-column before at last export time + self.data = self.data[self.data[:, 0].argsort()] + # Write data + np.savetxt( + self.filename, + self.data, + header=self.header, + delimiter=",", + comments="", + ) diff --git a/test/unit/test_exports/test_txt_export.py b/test/unit/test_exports/test_txt_export.py index 58f471081..41417cf76 100644 --- a/test/unit/test_exports/test_txt_export.py +++ b/test/unit/test_exports/test_txt_export.py @@ -36,14 +36,14 @@ def my_export(self, tmpdir): def test_file_exists(self, my_export, function): current_time = 1 my_export.function = function - my_export.write(current_time=current_time, steady=False) + my_export.write(current_time=current_time, final_time=None) assert os.path.exists(my_export.filename) def test_file_doesnt_exist(self, my_export, function): current_time = 10 my_export.function = function - my_export.write(current_time=current_time, steady=False) + my_export.write(current_time=current_time, final_time=None) assert not os.path.exists(my_export.filename) @@ -57,14 +57,14 @@ def test_create_folder(self, my_export, function): + "/folder2" + my_export.filename[slash_indx:] ) - my_export.write(current_time=current_time, steady=False) + my_export.write(current_time=current_time, final_time=None) assert os.path.exists(my_export.filename) def test_subspace(self, my_export, function_subspace): current_time = 1 my_export.function = function_subspace - my_export.write(current_time=current_time, steady=False) + my_export.write(current_time=current_time, final_time=None) assert os.path.exists(my_export.filename) @@ -99,26 +99,3 @@ def test_false(self, my_export): assert not my_export.is_it_time_to_export(0) assert not my_export.is_it_time_to_export(5) assert not my_export.is_it_time_to_export(1.5) - - -class TestWhenIsNextTime: - @pytest.fixture - def my_export(self, tmpdir): - d = tmpdir.mkdir("test_folder") - my_export = TXTExport( - "solute", - times=[1, 2, 3], - filename="{}/solute_label.txt".format(str(Path(d))), - ) - - return my_export - - def test_there_is_a_next_time(self, my_export): - assert my_export.when_is_next_time(2) == 3 - assert my_export.when_is_next_time(1) == 2 - assert my_export.when_is_next_time(0) == 1 - assert my_export.when_is_next_time(0.5) == 1 - - def test_last(self, my_export): - assert my_export.when_is_next_time(3) is None - assert my_export.when_is_next_time(4) is None diff --git a/test/unit/test_exports/test_txt_exports.py b/test/unit/test_exports/test_txt_exports.py deleted file mode 100644 index 9d0317f50..000000000 --- a/test/unit/test_exports/test_txt_exports.py +++ /dev/null @@ -1,28 +0,0 @@ -from festim import TXTExports -import pytest -from pathlib import Path - - -class TestWrite: - @pytest.fixture(params=[[1, 2, 3], None]) - def my_export(self, tmpdir, request): - d = tmpdir.mkdir("test_folder") - my_export = TXTExports( - fields=["solute", "T"], - filenames=[ - "{}/solute_label.txt".format(str(Path(d))), - "{}/T_label.txt".format(str(Path(d))), - ], - times=request.param, - ) - - return my_export - - def test_txt_exports_times(self, my_export): - for export in my_export.exports: - assert export.times == my_export.times - - -def test_error_when_fields_and_labels_have_different_lengths(): - with pytest.raises(ValueError, match="Number of fields to be exported"): - TXTExports(["solute", "T"], ["solute_label.txt"], [1]) From 35dc3a7a140cb56c88cf63dab258f4b3e8c7bb4a Mon Sep 17 00:00:00 2001 From: Vladimir Kulagin Date: Fri, 20 Sep 2024 00:23:51 +0300 Subject: [PATCH 2/8] filter duplicates depending on chemical_pot flag --- festim/exports/exports.py | 4 +- festim/exports/txt_export.py | 52 +++++++++++++--- festim/generic_simulation.py | 7 ++- test/unit/test_exports/test_txt_export.py | 72 +++++++++++++++++++++-- 4 files changed, 120 insertions(+), 15 deletions(-) diff --git a/festim/exports/exports.py b/festim/exports/exports.py index c3d61480a..915f44199 100644 --- a/festim/exports/exports.py +++ b/festim/exports/exports.py @@ -71,7 +71,7 @@ def _validate_export(self, value): return value raise TypeError("festim.Exports must be a list of festim.Export") - def write(self, label_to_function, dx): + def write(self, label_to_function, dx, materials, chemical_pot): """writes to file Args: @@ -129,7 +129,7 @@ def write(self, label_to_function, dx): label_to_function[export.field], self.V_DG1 ) export.function = label_to_function[export.field] - export.write(self.t, self.final_time) + export.write(self.t, self.final_time, materials, chemical_pot) self.nb_iterations += 1 def initialise_derived_quantities(self, dx, ds, materials): diff --git a/festim/exports/txt_export.py b/festim/exports/txt_export.py index b306b31e3..9144c2b97 100644 --- a/festim/exports/txt_export.py +++ b/festim/exports/txt_export.py @@ -71,11 +71,44 @@ def is_last(self, current_time, final_time): return True return False - def write(self, current_time, final_time): - # create a DG1 functionspace - V_DG1 = f.FunctionSpace(self.function.function_space().mesh(), "DG", 1) + def filter_duplicates(self, data, materials): + x = data[:, 0] - solution = f.project(self.function, V_DG1) + # Collect all borders + borders = [] + for material in materials: + for border in material.borders: + borders.append(border) + borders = np.unique(borders) + + # Find indices of the closest duplicates to interfaces + border_indx = [] + for border in borders: + closest_indx = np.abs(x - border).argmin() + closest_x = x[closest_indx] + for ind in np.where(x == closest_x)[0]: + border_indx.append(ind) + + # Find indices of first elements in duplicated pairs and mesh borders + _, unique_indx = np.unique(x, return_index=True) + + # Combine both arrays of indices + combined_indx = np.concatenate([border_indx, unique_indx]) + + # Sort unique indices to return a slice + combined_indx = sorted(np.unique(combined_indx)) + + return data[combined_indx, :] + + def write(self, current_time, final_time, materials, chemical_pot): + # create a DG1 functionspace if chemical_pot is True + # else create a CG1 functionspace + if chemical_pot: + V = f.FunctionSpace(self.function.function_space().mesh(), "DG", 1) + else: + V = f.FunctionSpace(self.function.function_space().mesh(), "CG", 1) + + solution = f.project(self.function, V) solution_column = np.transpose(solution.vector()[:]) if self.is_it_time_to_export(current_time): # if the directory doesn't exist @@ -84,7 +117,7 @@ def write(self, current_time, final_time): if not os.path.exists(dirname): os.makedirs(dirname, exist_ok=True) - # write data if steady or it is the first time to export + # create header if steady or it is the first time to export # else append new column to the existing file if final_time is None or self._first_time: if final_time is None: @@ -92,7 +125,7 @@ def write(self, current_time, final_time): else: self.header = f"x,t={format(current_time, self.header_format)}s" - x = f.interpolate(f.Expression("x[0]", degree=1), V_DG1) + x = f.interpolate(f.Expression("x[0]", degree=1), V) x_column = np.transpose([x.vector()[:]]) self.data = np.column_stack([x_column, solution_column]) self._first_time = False @@ -106,8 +139,13 @@ def write(self, current_time, final_time): self.write_at_last and self.is_last(current_time, final_time) ) or not self.write_at_last: if self.is_last(current_time, final_time): - # Sort data by the x-column before at last export time + # Sort data by the x-column before the last export time self.data = self.data[self.data[:, 0].argsort()] + + # Filter duplicates if chemical_pot is True + if chemical_pot: + self.data = self.filter_duplicates(self.data, materials) + # Write data np.savetxt( self.filename, diff --git a/festim/generic_simulation.py b/festim/generic_simulation.py index f8766d1ee..96fa71148 100644 --- a/festim/generic_simulation.py +++ b/festim/generic_simulation.py @@ -503,7 +503,12 @@ def run_post_processing(self): self.update_post_processing_solutions() self.exports.t = self.t - self.exports.write(self.label_to_function, self.mesh.dx) + self.exports.write( + self.label_to_function, + self.mesh.dx, + self.materials, + self.settings.chemical_pot, + ) def update_post_processing_solutions(self): """Creates the post-processing functions by splitting self.u. Projects diff --git a/test/unit/test_exports/test_txt_export.py b/test/unit/test_exports/test_txt_export.py index 41417cf76..134d88255 100644 --- a/test/unit/test_exports/test_txt_export.py +++ b/test/unit/test_exports/test_txt_export.py @@ -1,7 +1,8 @@ -from festim import TXTExport, Stepsize +from festim import TXTExport, Stepsize, Material import fenics as f import os import pytest +import numpy as np from pathlib import Path @@ -36,14 +37,24 @@ def my_export(self, tmpdir): def test_file_exists(self, my_export, function): current_time = 1 my_export.function = function - my_export.write(current_time=current_time, final_time=None) + my_export.write( + current_time=current_time, + final_time=None, + materials=None, + chemical_pot=False, + ) assert os.path.exists(my_export.filename) def test_file_doesnt_exist(self, my_export, function): current_time = 10 my_export.function = function - my_export.write(current_time=current_time, final_time=None) + my_export.write( + current_time=current_time, + final_time=None, + materials=None, + chemical_pot=False, + ) assert not os.path.exists(my_export.filename) @@ -57,14 +68,24 @@ def test_create_folder(self, my_export, function): + "/folder2" + my_export.filename[slash_indx:] ) - my_export.write(current_time=current_time, final_time=None) + my_export.write( + current_time=current_time, + final_time=None, + materials=None, + chemical_pot=False, + ) assert os.path.exists(my_export.filename) def test_subspace(self, my_export, function_subspace): current_time = 1 my_export.function = function_subspace - my_export.write(current_time=current_time, final_time=None) + my_export.write( + current_time=current_time, + final_time=None, + materials=None, + chemical_pot=False, + ) assert os.path.exists(my_export.filename) @@ -76,6 +97,47 @@ def test_error_filename_not_a_str(self, my_export): with pytest.raises(TypeError, match="filename must be a string"): my_export.filename = 2 + def test_sorted_by_x(self, my_export, function): + """Checks that the exported data is sorted by x""" + current_time = 1 + my_export.function = function + my_export.write( + current_time=current_time, + final_time=None, + materials=None, + chemical_pot=False, + ) + assert (np.diff(my_export.data[:, 0]) >= 0).all() + + @pytest.mark.parametrize( + "materials,chemical_pot,export_len", + [ + (None, False, 11), + ( + [ + Material(id=1, D_0=1, E_D=0, S_0=1, E_S=0, borders=[0, 0.5]), + Material(id=2, D_0=2, E_D=0, S_0=2, E_S=0, borders=[0.5, 1]), + ], + True, + 12, # + 1 duplicate near the interface + ), + ], + ) + def test_duplicates(self, materials, chemical_pot, export_len, my_export, function): + """ + Checks that the exported data does not contain duplicates + except those near interfaces + """ + current_time = 1 + my_export.function = function + my_export.write( + current_time=current_time, + final_time=None, + materials=materials, + chemical_pot=chemical_pot, + ) + assert len(my_export.data) == export_len + class TestIsItTimeToExport: @pytest.fixture From 72648b20ec4a795334adcab3f751b08a7af24ede Mon Sep 17 00:00:00 2001 From: KulaginVladimir Date: Mon, 7 Oct 2024 02:40:08 +0300 Subject: [PATCH 3/8] refactor filtering --- festim/exports/exports.py | 4 +- festim/exports/txt_export.py | 125 ++++++++++++---------- festim/generic_simulation.py | 58 +++++----- test/unit/test_exports/test_txt_export.py | 42 ++++---- 4 files changed, 125 insertions(+), 104 deletions(-) diff --git a/festim/exports/exports.py b/festim/exports/exports.py index 915f44199..c3d61480a 100644 --- a/festim/exports/exports.py +++ b/festim/exports/exports.py @@ -71,7 +71,7 @@ def _validate_export(self, value): return value raise TypeError("festim.Exports must be a list of festim.Export") - def write(self, label_to_function, dx, materials, chemical_pot): + def write(self, label_to_function, dx): """writes to file Args: @@ -129,7 +129,7 @@ def write(self, label_to_function, dx, materials, chemical_pot): label_to_function[export.field], self.V_DG1 ) export.function = label_to_function[export.field] - export.write(self.t, self.final_time, materials, chemical_pot) + export.write(self.t, self.final_time) self.nb_iterations += 1 def initialise_derived_quantities(self, dx, ds, materials): diff --git a/festim/exports/txt_export.py b/festim/exports/txt_export.py index 9144c2b97..7884c69b4 100644 --- a/festim/exports/txt_export.py +++ b/festim/exports/txt_export.py @@ -13,11 +13,20 @@ class TXTExport(festim.Export): field (str): the exported field ("solute", "1", "retention", "T"...) filename (str): the filename (must end with .txt). + write_at_last (bool): if True, the data will be exported at + the last export time. Otherwise, the data will be exported + at each export time. Defaults to False. times (list, optional): if provided, the field will be exported at these timesteps. Otherwise exports at all timesteps. Defaults to None. header_format (str, optional): the format of column headers. Defautls to ".2e". + + Attributes: + data (np.array): the data array of the exported field. The first column + is the mesh vertices. Each next column is the field profile at the specific + export time. + header (str): the header of the exported file. """ def __init__( @@ -31,10 +40,11 @@ def __init__( self.filename = filename self.write_at_last = write_at_last self.header_format = header_format - self._first_time = True self.data = None self.header = None + self._unique_indices = None + self._V = None @property def filename(self): @@ -71,80 +81,81 @@ def is_last(self, current_time, final_time): return True return False - def filter_duplicates(self, data, materials): - x = data[:, 0] - - # Collect all borders - borders = [] - for material in materials: - for border in material.borders: - borders.append(border) - borders = np.unique(borders) - - # Find indices of the closest duplicates to interfaces - border_indx = [] - for border in borders: - closest_indx = np.abs(x - border).argmin() - closest_x = x[closest_indx] - for ind in np.where(x == closest_x)[0]: - border_indx.append(ind) + def initialise_TXTExport(self, mesh, project_to_DG=False, materials=None): - # Find indices of first elements in duplicated pairs and mesh borders - _, unique_indx = np.unique(x, return_index=True) - - # Combine both arrays of indices - combined_indx = np.concatenate([border_indx, unique_indx]) + if project_to_DG: + self._V = f.FunctionSpace(mesh, "DG", 1) + else: + self._V = f.FunctionSpace(mesh, "CG", 1) + + x = f.interpolate(f.Expression("x[0]", degree=1), self._V) + x_column = np.transpose([x.vector()[:]]) + + # if chemical_pot is True or trap_element_type is DG, get indices of duplicates near interfaces + # and indices of the first elements from a pair of duplicates otherwise + if project_to_DG: + # Collect all borders + borders = [] + for material in materials: + if material.borders: + for border in material.borders: + borders.append(border) + borders = np.unique(borders) + + # Find indices of the closest duplicates to interfaces + border_indices = [] + for border in borders: + closest_indx = np.abs(x_column - border).argmin() + closest_x = x_column[closest_indx] + for ind in np.where(x_column == closest_x)[0]: + border_indices.append(ind) + + # Find indices of first elements in duplicated pairs and mesh borders + _, mesh_indices = np.unique(x_column, return_index=True) + + # Get unique indices from both arrays preserving the order in unsorted x-array + unique_indices = [] + for indx in np.argsort(x_column, axis=0)[:, 0]: + if (indx in mesh_indices) or (indx in border_indices): + unique_indices.append(indx) + + self._unique_indices = np.array(unique_indices) - # Sort unique indices to return a slice - combined_indx = sorted(np.unique(combined_indx)) + else: + # Get list of unique indices as integers + self._unique_indices = np.argsort(x_column, axis=0)[:, 0] - return data[combined_indx, :] + self.data = x_column[self._unique_indices] + self.header = "x" - def write(self, current_time, final_time, materials, chemical_pot): - # create a DG1 functionspace if chemical_pot is True - # else create a CG1 functionspace - if chemical_pot: - V = f.FunctionSpace(self.function.function_space().mesh(), "DG", 1) - else: - V = f.FunctionSpace(self.function.function_space().mesh(), "CG", 1) + def write(self, current_time, final_time): - solution = f.project(self.function, V) - solution_column = np.transpose(solution.vector()[:]) if self.is_it_time_to_export(current_time): + solution = f.project(self.function, self._V) + solution_column = np.transpose(solution.vector()[:]) + # if the directory doesn't exist # create it dirname = os.path.dirname(self.filename) if not os.path.exists(dirname): os.makedirs(dirname, exist_ok=True) - # create header if steady or it is the first time to export - # else append new column to the existing file - if final_time is None or self._first_time: - if final_time is None: - self.header = "x,t=steady" - else: - self.header = f"x,t={format(current_time, self.header_format)}s" - - x = f.interpolate(f.Expression("x[0]", degree=1), V) - x_column = np.transpose([x.vector()[:]]) - self.data = np.column_stack([x_column, solution_column]) - self._first_time = False + # if steady, add the corresponding label + # else append new export time to the header + steady = final_time is None + if steady: + self.header += ",t=steady" else: - # Update the header self.header += f",t={format(current_time, self.header_format)}s" - # Add new column - self.data = np.column_stack([self.data, solution_column]) + + # Add new column of filtered and sorted data + self.data = np.column_stack( + [self.data, solution_column[self._unique_indices]] + ) if ( self.write_at_last and self.is_last(current_time, final_time) ) or not self.write_at_last: - if self.is_last(current_time, final_time): - # Sort data by the x-column before the last export time - self.data = self.data[self.data[:, 0].argsort()] - - # Filter duplicates if chemical_pot is True - if chemical_pot: - self.data = self.filter_duplicates(self.data, materials) # Write data np.savetxt( diff --git a/festim/generic_simulation.py b/festim/generic_simulation.py index 96fa71148..2db764373 100644 --- a/festim/generic_simulation.py +++ b/festim/generic_simulation.py @@ -357,11 +357,11 @@ def initialise(self): self.h_transport_problem.initialise(self.mesh, self.materials, self.dt) - # raise warning if the derived quantities don't match the type of mesh - # eg. SurfaceFlux is used with cylindrical mesh for export in self.exports: if isinstance(export, festim.DerivedQuantities): for q in export: + # raise warning if the derived quantities don't match the type of mesh + # eg. SurfaceFlux is used with cylindrical mesh if self.mesh.type not in q.allowed_meshes: warnings.warn( f"{type(q)} may not work as intended for {self.mesh.type} meshes" @@ -381,31 +381,41 @@ def initialise(self): f"SurfaceKinetics boundary condition must be defined on surface {q.surface} to export data with festim.AdsorbedHydrogen" ) - self.exports.initialise_derived_quantities( - self.mesh.dx, self.mesh.ds, self.materials - ) - - # needed to ensure that data is actually exported at TXTExport.times - # see issue 675 - for export in self.exports: - if isinstance(export, festim.TXTExport) and export.times: - if not self.dt.milestones: - self.dt.milestones = [] - for time in export.times: - if time not in self.dt.milestones: - msg = "To ensure that TXTExport exports data at the desired times " - msg += "TXTExport.times are added to milestones" - warnings.warn(msg) - self.dt.milestones.append(time) - self.dt.milestones.sort() - - # set Soret to True for SurfaceFlux quantities - if isinstance(export, festim.DerivedQuantities): - for q in export: + # set Soret to True for SurfaceFlux quantities if isinstance(q, festim.SurfaceFlux): q.soret = self.settings.soret q.T = self.T.T + if isinstance(export, festim.TXTExport): + # pre-process data depending on the chemical potential flag, trap element type, + # and material borders + project_to_DG = ( + self.settings.chemical_pot + or self.settings.traps_element_type == "DG" + ) + export.initialise_TXTExport( + self.mesh.mesh, + project_to_DG, + self.materials, + ) + + # needed to ensure that data is actually exported at TXTExport.times + # see issue 675 + if export.times: + if not self.dt.milestones: + self.dt.milestones = [] + for time in export.times: + if time not in self.dt.milestones: + msg = "To ensure that TXTExport exports data at the desired times " + msg += "TXTExport.times are added to milestones" + warnings.warn(msg) + self.dt.milestones.append(time) + self.dt.milestones.sort() + + self.exports.initialise_derived_quantities( + self.mesh.dx, self.mesh.ds, self.materials + ) + def run(self, completion_tone=False): """Runs the model. @@ -506,8 +516,6 @@ def run_post_processing(self): self.exports.write( self.label_to_function, self.mesh.dx, - self.materials, - self.settings.chemical_pot, ) def update_post_processing_solutions(self): diff --git a/test/unit/test_exports/test_txt_export.py b/test/unit/test_exports/test_txt_export.py index 134d88255..5210bf3f3 100644 --- a/test/unit/test_exports/test_txt_export.py +++ b/test/unit/test_exports/test_txt_export.py @@ -1,4 +1,4 @@ -from festim import TXTExport, Stepsize, Material +from festim import TXTExport, Material import fenics as f import os import pytest @@ -7,6 +7,12 @@ class TestWrite: + @pytest.fixture + def mesh(self): + mesh = f.UnitIntervalMesh(10) + + return mesh + @pytest.fixture def function(self): mesh = f.UnitIntervalMesh(10) @@ -34,31 +40,29 @@ def my_export(self, tmpdir): return my_export - def test_file_exists(self, my_export, function): + def test_file_exists(self, my_export, function, mesh): current_time = 1 my_export.function = function + my_export.initialise_TXTExport(mesh) my_export.write( current_time=current_time, final_time=None, - materials=None, - chemical_pot=False, ) assert os.path.exists(my_export.filename) - def test_file_doesnt_exist(self, my_export, function): + def test_file_doesnt_exist(self, my_export, function, mesh): current_time = 10 my_export.function = function + my_export.initialise_TXTExport(mesh) my_export.write( current_time=current_time, final_time=None, - materials=None, - chemical_pot=False, ) assert not os.path.exists(my_export.filename) - def test_create_folder(self, my_export, function): + def test_create_folder(self, my_export, function, mesh): """Checks that write() creates the folder if it doesn't exist""" current_time = 1 my_export.function = function @@ -68,23 +72,21 @@ def test_create_folder(self, my_export, function): + "/folder2" + my_export.filename[slash_indx:] ) + my_export.initialise_TXTExport(mesh) my_export.write( current_time=current_time, final_time=None, - materials=None, - chemical_pot=False, ) assert os.path.exists(my_export.filename) - def test_subspace(self, my_export, function_subspace): + def test_subspace(self, my_export, function_subspace, mesh): current_time = 1 my_export.function = function_subspace + my_export.initialise_TXTExport(mesh) my_export.write( current_time=current_time, final_time=None, - materials=None, - chemical_pot=False, ) assert os.path.exists(my_export.filename) @@ -97,20 +99,19 @@ def test_error_filename_not_a_str(self, my_export): with pytest.raises(TypeError, match="filename must be a string"): my_export.filename = 2 - def test_sorted_by_x(self, my_export, function): + def test_sorted_by_x(self, my_export, function, mesh): """Checks that the exported data is sorted by x""" current_time = 1 my_export.function = function + my_export.initialise_TXTExport(mesh) my_export.write( current_time=current_time, final_time=None, - materials=None, - chemical_pot=False, ) assert (np.diff(my_export.data[:, 0]) >= 0).all() @pytest.mark.parametrize( - "materials,chemical_pot,export_len", + "materials,project_to_DG,export_len", [ (None, False, 11), ( @@ -123,18 +124,19 @@ def test_sorted_by_x(self, my_export, function): ), ], ) - def test_duplicates(self, materials, chemical_pot, export_len, my_export, function): + def test_duplicates( + self, materials, project_to_DG, export_len, my_export, function, mesh + ): """ Checks that the exported data does not contain duplicates except those near interfaces """ current_time = 1 my_export.function = function + my_export.initialise_TXTExport(mesh, project_to_DG, materials) my_export.write( current_time=current_time, final_time=None, - materials=materials, - chemical_pot=chemical_pot, ) assert len(my_export.data) == export_len From 7ca79976e5bb685d9a69a652179664fc39eaeeba Mon Sep 17 00:00:00 2001 From: KulaginVladimir Date: Mon, 7 Oct 2024 03:17:53 +0300 Subject: [PATCH 4/8] test is_last --- test/unit/test_exports/test_txt_export.py | 26 +++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/test/unit/test_exports/test_txt_export.py b/test/unit/test_exports/test_txt_export.py index 5210bf3f3..177b6c853 100644 --- a/test/unit/test_exports/test_txt_export.py +++ b/test/unit/test_exports/test_txt_export.py @@ -163,3 +163,29 @@ def test_false(self, my_export): assert not my_export.is_it_time_to_export(0) assert not my_export.is_it_time_to_export(5) assert not my_export.is_it_time_to_export(1.5) + + +class TestIsLast: + @pytest.fixture + def my_export(self, tmpdir): + d = tmpdir.mkdir("test_folder") + my_export = TXTExport( + "solute", + filename="{}/solute_label.txt".format(str(Path(d))), + ) + + return my_export + + def test_final_time_none(self, my_export): + assert my_export.is_last(1, None) == True + + @pytest.mark.parametrize("current_time,output", [(1, False), (3, False), (5, True)]) + def test_times_not_none(self, current_time, output, my_export): + final_time = 5 + assert my_export.is_last(current_time, final_time) == output + + @pytest.mark.parametrize("current_time,output", [(1, False), (2, False), (3, True)]) + def test_times_not_none(self, current_time, output, my_export): + my_export.times = [1, 2, 3] + final_time = 5 + assert my_export.is_last(current_time, final_time) == output From 4903d7edbef943be31e1e3facd8e044058ac424f Mon Sep 17 00:00:00 2001 From: KulaginVladimir Date: Mon, 7 Oct 2024 03:37:29 +0300 Subject: [PATCH 5/8] fix typo --- test/unit/test_exports/test_txt_export.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/unit/test_exports/test_txt_export.py b/test/unit/test_exports/test_txt_export.py index 177b6c853..e34b08e7b 100644 --- a/test/unit/test_exports/test_txt_export.py +++ b/test/unit/test_exports/test_txt_export.py @@ -176,11 +176,11 @@ def my_export(self, tmpdir): return my_export - def test_final_time_none(self, my_export): + def test_final_time_is_none(self, my_export): assert my_export.is_last(1, None) == True @pytest.mark.parametrize("current_time,output", [(1, False), (3, False), (5, True)]) - def test_times_not_none(self, current_time, output, my_export): + def test_times_is_none(self, current_time, output, my_export): final_time = 5 assert my_export.is_last(current_time, final_time) == output From 9ef34ad938a1446d4fcc3e89488893932636b605 Mon Sep 17 00:00:00 2001 From: KulaginVladimir Date: Mon, 7 Oct 2024 22:43:34 +0300 Subject: [PATCH 6/8] added docstrings --- docs/source/api/exports.rst | 4 -- festim/exports/txt_export.py | 67 ++++++++++++++++++++--- festim/generic_simulation.py | 2 +- test/unit/test_exports/test_txt_export.py | 12 ++-- 4 files changed, 66 insertions(+), 19 deletions(-) diff --git a/docs/source/api/exports.rst b/docs/source/api/exports.rst index abde0afbc..0bf38fdfc 100644 --- a/docs/source/api/exports.rst +++ b/docs/source/api/exports.rst @@ -15,10 +15,6 @@ Exports :members: :show-inheritance: -.. autoclass:: TXTExports - :members: - :show-inheritance: - .. autoclass:: TrapDensityXDMF :members: :show-inheritance: diff --git a/festim/exports/txt_export.py b/festim/exports/txt_export.py index 7884c69b4..0dadd4614 100644 --- a/festim/exports/txt_export.py +++ b/festim/exports/txt_export.py @@ -27,6 +27,8 @@ class TXTExport(festim.Export): is the mesh vertices. Each next column is the field profile at the specific export time. header (str): the header of the exported file. + V (fenics.FunctionSpace): the vector-function space for the exported field. + """ def __init__( @@ -43,8 +45,8 @@ def __init__( self.data = None self.header = None + self.V = None self._unique_indices = None - self._V = None @property def filename(self): @@ -60,6 +62,17 @@ def filename(self, value): self._filename = value def is_it_time_to_export(self, current_time): + """ + Checks if the exported field should be written to a file or not + based on the current time and the TXTExport.times + + Args: + current_time (float): the current simulation time + + Returns: + bool: True if the exported field should be written to a file, else False + """ + if self.times is None: return True for time in self.times: @@ -68,6 +81,21 @@ def is_it_time_to_export(self, current_time): return False def is_last(self, current_time, final_time): + """ + Checks if the current simulation step equals to the last export time. + based on the final simulation time, TXTExport.times, and the current time + + Args: + current_time (float): the current simulation time. + final_time (float, None): the final simulation time. + + Returns: + bool: True if simulation is steady (final_time is None), if TXTExport.times + are not provided and the current time equals to the final time, or if + TXTExport.times are provided and the current time equals to the last time in + TXTExport.times, else False. + """ + if final_time is None: # write if steady return True @@ -81,17 +109,32 @@ def is_last(self, current_time, final_time): return True return False - def initialise_TXTExport(self, mesh, project_to_DG=False, materials=None): + def initialise(self, mesh, project_to_DG=False, materials=None): + """ + Initialises TXTExport. Depending on the project_to_DG flag, defines a function space (DG1 or CG1) + for projection of the exported field. After that, an unsorted array of mesh vertices is created for export. + The array is then used to obtain indices of sorted elements for the data export. + + .. note:: + If DG1 is used, the duplicated vertices in the array are filtered except those near interfaces, + The interfaces are defined by ``material.borders`` in the ``Materials`` list. + + Args: + mesh (fenics.Mesh): the mesh. + project_to_DG (bool): if True, the exported field is projected to a DG1 function space. + Defaults to False. + materials (festim.Materials): the materials. Defaults to None. + """ if project_to_DG: - self._V = f.FunctionSpace(mesh, "DG", 1) + self.V = f.FunctionSpace(mesh, "DG", 1) else: - self._V = f.FunctionSpace(mesh, "CG", 1) + self.V = f.FunctionSpace(mesh, "CG", 1) - x = f.interpolate(f.Expression("x[0]", degree=1), self._V) + x = f.interpolate(f.Expression("x[0]", degree=1), self.V) x_column = np.transpose([x.vector()[:]]) - # if chemical_pot is True or trap_element_type is DG, get indices of duplicates near interfaces + # if project_to_DG is True, get indices of duplicates near interfaces # and indices of the first elements from a pair of duplicates otherwise if project_to_DG: # Collect all borders @@ -122,16 +165,24 @@ def initialise_TXTExport(self, mesh, project_to_DG=False, materials=None): self._unique_indices = np.array(unique_indices) else: - # Get list of unique indices as integers + # Get list of unique indices self._unique_indices = np.argsort(x_column, axis=0)[:, 0] self.data = x_column[self._unique_indices] self.header = "x" def write(self, current_time, final_time): + """ + Modifies the header and writes the data to a file depending on + the current and the final times of a simulation. + + Args: + current_time (float): the current simulation time. + final_time (float, None): the final simulation time. + """ if self.is_it_time_to_export(current_time): - solution = f.project(self.function, self._V) + solution = f.project(self.function, self.V) solution_column = np.transpose(solution.vector()[:]) # if the directory doesn't exist diff --git a/festim/generic_simulation.py b/festim/generic_simulation.py index 2db764373..d8167d955 100644 --- a/festim/generic_simulation.py +++ b/festim/generic_simulation.py @@ -393,7 +393,7 @@ def initialise(self): self.settings.chemical_pot or self.settings.traps_element_type == "DG" ) - export.initialise_TXTExport( + export.initialise( self.mesh.mesh, project_to_DG, self.materials, diff --git a/test/unit/test_exports/test_txt_export.py b/test/unit/test_exports/test_txt_export.py index e34b08e7b..fe9493e78 100644 --- a/test/unit/test_exports/test_txt_export.py +++ b/test/unit/test_exports/test_txt_export.py @@ -43,7 +43,7 @@ def my_export(self, tmpdir): def test_file_exists(self, my_export, function, mesh): current_time = 1 my_export.function = function - my_export.initialise_TXTExport(mesh) + my_export.initialise(mesh) my_export.write( current_time=current_time, final_time=None, @@ -54,7 +54,7 @@ def test_file_exists(self, my_export, function, mesh): def test_file_doesnt_exist(self, my_export, function, mesh): current_time = 10 my_export.function = function - my_export.initialise_TXTExport(mesh) + my_export.initialise(mesh) my_export.write( current_time=current_time, final_time=None, @@ -72,7 +72,7 @@ def test_create_folder(self, my_export, function, mesh): + "/folder2" + my_export.filename[slash_indx:] ) - my_export.initialise_TXTExport(mesh) + my_export.initialise(mesh) my_export.write( current_time=current_time, final_time=None, @@ -83,7 +83,7 @@ def test_create_folder(self, my_export, function, mesh): def test_subspace(self, my_export, function_subspace, mesh): current_time = 1 my_export.function = function_subspace - my_export.initialise_TXTExport(mesh) + my_export.initialise(mesh) my_export.write( current_time=current_time, final_time=None, @@ -103,7 +103,7 @@ def test_sorted_by_x(self, my_export, function, mesh): """Checks that the exported data is sorted by x""" current_time = 1 my_export.function = function - my_export.initialise_TXTExport(mesh) + my_export.initialise(mesh) my_export.write( current_time=current_time, final_time=None, @@ -133,7 +133,7 @@ def test_duplicates( """ current_time = 1 my_export.function = function - my_export.initialise_TXTExport(mesh, project_to_DG, materials) + my_export.initialise(mesh, project_to_DG, materials) my_export.write( current_time=current_time, final_time=None, From f639d2a2ea1506d21c1f1fe65ba0dc55d7dd22d5 Mon Sep 17 00:00:00 2001 From: KulaginVladimir Date: Sun, 27 Oct 2024 23:19:43 +0300 Subject: [PATCH 7/8] added filter flag --- festim/exports/txt_export.py | 30 +++++++++++++++++------ test/unit/test_exports/test_txt_export.py | 18 +++++++++++--- 2 files changed, 37 insertions(+), 11 deletions(-) diff --git a/festim/exports/txt_export.py b/festim/exports/txt_export.py index 0dadd4614..f148716f6 100644 --- a/festim/exports/txt_export.py +++ b/festim/exports/txt_export.py @@ -13,12 +13,15 @@ class TXTExport(festim.Export): field (str): the exported field ("solute", "1", "retention", "T"...) filename (str): the filename (must end with .txt). - write_at_last (bool): if True, the data will be exported at - the last export time. Otherwise, the data will be exported - at each export time. Defaults to False. times (list, optional): if provided, the field will be exported at these timesteps. Otherwise exports at all timesteps. Defaults to None. + filter (bool): if True and the field is projected to a DG function space, + the duplicated vertices in the output file array are filtered except those near interfaces. + Defaults to True. + write_at_last (bool): if True, the data will be exported at + the last export time. Otherwise, the data will be exported + at each export time. Defaults to False. header_format (str, optional): the format of column headers. Defautls to ".2e". @@ -29,10 +32,20 @@ class TXTExport(festim.Export): header (str): the header of the exported file. V (fenics.FunctionSpace): the vector-function space for the exported field. + .. note:: + The exported field is projected to DG if conservation of chemical potential is considered or + ``traps_element_type`` is "DG". + """ def __init__( - self, field, filename, times=None, write_at_last=False, header_format=".2e" + self, + field, + filename, + times=None, + filter=True, + write_at_last=False, + header_format=".2e", ) -> None: super().__init__(field=field) if times: @@ -40,6 +53,7 @@ def __init__( else: self.times = times self.filename = filename + self.filter = filter self.write_at_last = write_at_last self.header_format = header_format @@ -116,7 +130,7 @@ def initialise(self, mesh, project_to_DG=False, materials=None): The array is then used to obtain indices of sorted elements for the data export. .. note:: - If DG1 is used, the duplicated vertices in the array are filtered except those near interfaces, + If DG1 is used and filter flag is True, the duplicated vertices in the array are filtered except those near interfaces, The interfaces are defined by ``material.borders`` in the ``Materials`` list. Args: @@ -134,9 +148,9 @@ def initialise(self, mesh, project_to_DG=False, materials=None): x = f.interpolate(f.Expression("x[0]", degree=1), self.V) x_column = np.transpose([x.vector()[:]]) - # if project_to_DG is True, get indices of duplicates near interfaces + # if filter is True, get indices of duplicates near interfaces # and indices of the first elements from a pair of duplicates otherwise - if project_to_DG: + if project_to_DG and self.filter: # Collect all borders borders = [] for material in materials: @@ -165,7 +179,7 @@ def initialise(self, mesh, project_to_DG=False, materials=None): self._unique_indices = np.array(unique_indices) else: - # Get list of unique indices + # Get list of sorted indices self._unique_indices = np.argsort(x_column, axis=0)[:, 0] self.data = x_column[self._unique_indices] diff --git a/test/unit/test_exports/test_txt_export.py b/test/unit/test_exports/test_txt_export.py index fe9493e78..1d86f32c0 100644 --- a/test/unit/test_exports/test_txt_export.py +++ b/test/unit/test_exports/test_txt_export.py @@ -111,21 +111,31 @@ def test_sorted_by_x(self, my_export, function, mesh): assert (np.diff(my_export.data[:, 0]) >= 0).all() @pytest.mark.parametrize( - "materials,project_to_DG,export_len", + "materials,project_to_DG,filter,export_len", [ - (None, False, 11), + (None, False, False, 11), ( [ Material(id=1, D_0=1, E_D=0, S_0=1, E_S=0, borders=[0, 0.5]), Material(id=2, D_0=2, E_D=0, S_0=2, E_S=0, borders=[0.5, 1]), ], True, + True, 12, # + 1 duplicate near the interface ), + ( + [ + Material(id=1, D_0=1, E_D=0, S_0=1, E_S=0, borders=[0, 0.5]), + Material(id=2, D_0=2, E_D=0, S_0=2, E_S=0, borders=[0.5, 1]), + ], + True, + False, + 20, # 2 * (len_vertices - 1) + ), ], ) def test_duplicates( - self, materials, project_to_DG, export_len, my_export, function, mesh + self, materials, project_to_DG, filter, export_len, my_export, function, mesh ): """ Checks that the exported data does not contain duplicates @@ -133,11 +143,13 @@ def test_duplicates( """ current_time = 1 my_export.function = function + my_export.filter = filter my_export.initialise(mesh, project_to_DG, materials) my_export.write( current_time=current_time, final_time=None, ) + print(my_export._unique_indices) assert len(my_export.data) == export_len From 468611de0ef1ddb22621163281cab728c4249da5 Mon Sep 17 00:00:00 2001 From: KulaginVladimir Date: Sat, 21 Dec 2024 04:40:55 +0300 Subject: [PATCH 8/8] suggestions from review --- festim/exports/txt_export.py | 25 +++++++++++------------ test/unit/test_exports/test_txt_export.py | 2 +- 2 files changed, 13 insertions(+), 14 deletions(-) diff --git a/festim/exports/txt_export.py b/festim/exports/txt_export.py index f148716f6..6137b8b8f 100644 --- a/festim/exports/txt_export.py +++ b/festim/exports/txt_export.py @@ -78,7 +78,7 @@ def filename(self, value): def is_it_time_to_export(self, current_time): """ Checks if the exported field should be written to a file or not - based on the current time and the TXTExport.times + based on the current time and the ``TXTExport.times`` Args: current_time (float): the current simulation time @@ -97,17 +97,17 @@ def is_it_time_to_export(self, current_time): def is_last(self, current_time, final_time): """ Checks if the current simulation step equals to the last export time. - based on the final simulation time, TXTExport.times, and the current time + based on the final simulation time, ``TXTExport.times``, and the current time Args: current_time (float): the current simulation time. final_time (float, None): the final simulation time. Returns: - bool: True if simulation is steady (final_time is None), if TXTExport.times - are not provided and the current time equals to the final time, or if - TXTExport.times are provided and the current time equals to the last time in - TXTExport.times, else False. + bool: True if simulation is steady (final_time is None), if ``TXTExport.times`` are not + provided and the current time equals to the final time, or if + ``TXTExport.times`` are provided and the current time equals to the last time in + ``TXTExport.times``, else False. """ if final_time is None: @@ -115,22 +115,21 @@ def is_last(self, current_time, final_time): return True elif self.times is None: if np.isclose(current_time, final_time, atol=0): - # write at final time if exports at each timestep - return True - else: - if np.isclose(current_time, self.times[-1], atol=0): - # write at final time if exports at specific times + # write at the final time if exports at each timestep return True + elif np.isclose(current_time, self.times[-1], atol=0): + # write at the final time if exports at specific times + return True return False def initialise(self, mesh, project_to_DG=False, materials=None): """ - Initialises TXTExport. Depending on the project_to_DG flag, defines a function space (DG1 or CG1) + Initialises ``TXTExport``. Depending on the ``project_to_DG flag``, defines a function space (DG1 or CG1) for projection of the exported field. After that, an unsorted array of mesh vertices is created for export. The array is then used to obtain indices of sorted elements for the data export. .. note:: - If DG1 is used and filter flag is True, the duplicated vertices in the array are filtered except those near interfaces, + If DG1 is used and the ``filter`` flag is True, the duplicated vertices in the array are filtered except those near interfaces. The interfaces are defined by ``material.borders`` in the ``Materials`` list. Args: diff --git a/test/unit/test_exports/test_txt_export.py b/test/unit/test_exports/test_txt_export.py index 1d86f32c0..6b180b4f8 100644 --- a/test/unit/test_exports/test_txt_export.py +++ b/test/unit/test_exports/test_txt_export.py @@ -149,7 +149,7 @@ def test_duplicates( current_time=current_time, final_time=None, ) - print(my_export._unique_indices) + assert len(my_export.data) == export_len