From 35dc3a7a140cb56c88cf63dab258f4b3e8c7bb4a Mon Sep 17 00:00:00 2001 From: Vladimir Kulagin Date: Fri, 20 Sep 2024 00:23:51 +0300 Subject: [PATCH] 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