From 0e3004d0511aec56faa6356ba57e8b2598a6a2c9 Mon Sep 17 00:00:00 2001 From: KulaginVladimir Date: Mon, 9 Sep 2024 21:26:52 +0300 Subject: [PATCH 1/7] 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/7] 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/7] 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/7] 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/7] 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/7] 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/7] 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