From 97e036f29568118217d208a19db51d785931f60d Mon Sep 17 00:00:00 2001 From: kaelyndunnell Date: Fri, 15 Nov 2024 16:16:54 -0500 Subject: [PATCH] update fenicsx --- __init__.py | 117 ---- boundary_conditions/dirichlets/dc_imp.py | 95 --- boundary_conditions/dirichlets/sieverts_bc.py | 37 -- boundary_conditions/fluxes/convective_flux.py | 27 - .../fluxes/dissociation_flux.py | 33 -- boundary_conditions/fluxes/mass_flux.py | 30 - .../fluxes/surface_kinetics.py | 184 ------ concentration/mobile.py | 237 -------- concentration/theta.py | 122 ---- concentration/traps/extrinsic_trap.py | 161 ------ concentration/traps/trap.py | 211 ------- concentration/traps/traps.py | 133 ----- .../derived_quantities/adsorbed_hydrogen.py | 43 -- exports/derived_quantities/average_surface.py | 49 -- exports/derived_quantities/average_volume.py | 48 -- .../derived_quantities/derived_quantities.py | 302 ---------- .../derived_quantities/derived_quantity.py | 98 ---- exports/derived_quantities/hydrogen_flux.py | 26 - exports/derived_quantities/maximum_surface.py | 57 -- exports/derived_quantities/maximum_volume.py | 57 -- exports/derived_quantities/minimum_surface.py | 57 -- exports/derived_quantities/minimum_volume.py | 57 -- exports/derived_quantities/point_value.py | 46 -- exports/derived_quantities/surface_flux.py | 292 ---------- exports/derived_quantities/thermal_flux.py | 26 - exports/derived_quantities/total_surface.py | 60 -- exports/derived_quantities/total_volume.py | 60 -- exports/exports.py | 154 ----- exports/txt_export.py | 133 ----- exports/xdmf_export.py | 161 ------ generic_simulation.py | 542 ------------------ h_transport_problem.py | 388 ------------- initial_condition.py | 27 - materials/material.py | 83 --- materials/materials.py | 410 ------------- meshing/mesh.py | 37 -- meshing/mesh_1d.py | 79 --- meshing/mesh_from_vertices.py | 44 -- nonlinear_problem.py | 31 - settings.py | 76 --- sources/radioactive_decay.py | 35 -- stepsize.py | 139 ----- temperature/temperature.py | 56 -- temperature/temperature_from_xdmf.py | 52 -- temperature/temperature_solver.py | 275 --------- 45 files changed, 5387 deletions(-) delete mode 100644 __init__.py delete mode 100644 boundary_conditions/dirichlets/dc_imp.py delete mode 100644 boundary_conditions/dirichlets/sieverts_bc.py delete mode 100644 boundary_conditions/fluxes/convective_flux.py delete mode 100644 boundary_conditions/fluxes/dissociation_flux.py delete mode 100644 boundary_conditions/fluxes/mass_flux.py delete mode 100644 boundary_conditions/fluxes/surface_kinetics.py delete mode 100644 concentration/mobile.py delete mode 100644 concentration/theta.py delete mode 100644 concentration/traps/extrinsic_trap.py delete mode 100644 concentration/traps/trap.py delete mode 100644 concentration/traps/traps.py delete mode 100644 exports/derived_quantities/adsorbed_hydrogen.py delete mode 100644 exports/derived_quantities/average_surface.py delete mode 100644 exports/derived_quantities/average_volume.py delete mode 100644 exports/derived_quantities/derived_quantities.py delete mode 100644 exports/derived_quantities/derived_quantity.py delete mode 100644 exports/derived_quantities/hydrogen_flux.py delete mode 100644 exports/derived_quantities/maximum_surface.py delete mode 100644 exports/derived_quantities/maximum_volume.py delete mode 100644 exports/derived_quantities/minimum_surface.py delete mode 100644 exports/derived_quantities/minimum_volume.py delete mode 100644 exports/derived_quantities/point_value.py delete mode 100644 exports/derived_quantities/surface_flux.py delete mode 100644 exports/derived_quantities/thermal_flux.py delete mode 100644 exports/derived_quantities/total_surface.py delete mode 100644 exports/derived_quantities/total_volume.py delete mode 100644 exports/exports.py delete mode 100644 exports/txt_export.py delete mode 100644 exports/xdmf_export.py delete mode 100644 generic_simulation.py delete mode 100644 h_transport_problem.py delete mode 100644 initial_condition.py delete mode 100644 materials/material.py delete mode 100644 materials/materials.py delete mode 100644 meshing/mesh.py delete mode 100644 meshing/mesh_1d.py delete mode 100644 meshing/mesh_from_vertices.py delete mode 100644 nonlinear_problem.py delete mode 100644 settings.py delete mode 100644 sources/radioactive_decay.py delete mode 100644 stepsize.py delete mode 100644 temperature/temperature.py delete mode 100644 temperature/temperature_from_xdmf.py delete mode 100644 temperature/temperature_solver.py diff --git a/__init__.py b/__init__.py deleted file mode 100644 index 9a7bc7e59..000000000 --- a/__init__.py +++ /dev/null @@ -1,117 +0,0 @@ -try: - # Python 3.8+ - from importlib import metadata -except ImportError: - try: - import importlib_metadata as metadata - except ImportError: - __version__ = "unknown" - -try: - __version__ = metadata.version("FESTIM") -except Exception: - __version__ = "unknown" - - -import sympy as sp - -x, y, z, t = sp.symbols("x[0] x[1] x[2] t") -R = 8.314462618 # Gas constant J.mol-1.K-1 -k_B = 8.6173303e-5 # Boltzmann constant eV.K-1 - -from .helpers import ( - update_expressions, - kJmol_to_eV, - extract_xdmf_labels, - extract_xdmf_times, - as_constant, - as_expression, - as_constant_or_expression, -) - -from .meshing.mesh import Mesh -from .meshing.mesh_1d import Mesh1D -from .meshing.mesh_from_vertices import MeshFromVertices -from .meshing.mesh_from_xdmf import MeshFromXDMF - -from .temperature.temperature import Temperature -from .temperature.temperature_solver import HeatTransferProblem -from .temperature.temperature_from_xdmf import TemperatureFromXDMF - -from .boundary_conditions.boundary_condition import BoundaryCondition -from .boundary_conditions.dirichlets.dirichlet_bc import ( - DirichletBC, - BoundaryConditionTheta, - BoundaryConditionExpression, -) -from .boundary_conditions.dirichlets.dc_imp import ImplantationDirichlet -from .boundary_conditions.dirichlets.sieverts_bc import SievertsBC -from .boundary_conditions.dirichlets.henrys_bc import HenrysBC -from .boundary_conditions.dirichlets.custom_dc import CustomDirichlet - -from .boundary_conditions.fluxes.flux_bc import FluxBC -from .boundary_conditions.fluxes.recombination_flux import RecombinationFlux -from .boundary_conditions.fluxes.dissociation_flux import DissociationFlux -from .boundary_conditions.fluxes.convective_flux import ConvectiveFlux -from .boundary_conditions.fluxes.flux_custom import CustomFlux -from .boundary_conditions.fluxes.mass_flux import MassFlux -from .boundary_conditions.fluxes.surface_kinetics import SurfaceKinetics - -from .exports.exports import Exports -from .exports.export import Export -from .exports.xdmf_export import XDMFExport -from .exports.trap_density_xdmf import TrapDensityXDMF - -from .exports.derived_quantities.derived_quantity import ( - DerivedQuantity, - VolumeQuantity, - SurfaceQuantity, -) -from .exports.derived_quantities.surface_flux import ( - SurfaceFlux, - SurfaceFluxCylindrical, - SurfaceFluxSpherical, -) -from .exports.derived_quantities.hydrogen_flux import HydrogenFlux -from .exports.derived_quantities.thermal_flux import ThermalFlux -from .exports.derived_quantities.average_volume import AverageVolume -from .exports.derived_quantities.maximum_volume import MaximumVolume -from .exports.derived_quantities.minimum_volume import MinimumVolume -from .exports.derived_quantities.minimum_surface import MinimumSurface -from .exports.derived_quantities.maximum_surface import MaximumSurface -from .exports.derived_quantities.total_surface import TotalSurface -from .exports.derived_quantities.total_volume import TotalVolume -from .exports.derived_quantities.average_surface import AverageSurface -from .exports.derived_quantities.point_value import PointValue -from .exports.derived_quantities.adsorbed_hydrogen import AdsorbedHydrogen - -from .exports.derived_quantities.derived_quantities import DerivedQuantities - -from .exports.txt_export import TXTExport, TXTExports - - -from .settings import Settings -from .stepsize import Stepsize - -from .sources.source import Source -from .sources.source_implantation_flux import ImplantationFlux -from .sources.radioactive_decay import RadioactiveDecay - -from .materials.material import Material -from .materials.materials import Materials - -from .concentration.concentration import Concentration -from .initial_condition import InitialCondition -from .concentration.mobile import Mobile -from .nonlinear_problem import Problem -from .concentration.theta import Theta - -from .concentration.traps.trap import Trap -from .concentration.traps.traps import Traps -from .concentration.traps.extrinsic_trap import ExtrinsicTrapBase -from .concentration.traps.extrinsic_trap import ExtrinsicTrap -from .concentration.traps.neutron_induced_trap import NeutronInducedTrap - -from .h_transport_problem import HTransportProblem - -from .generic_simulation import Simulation diff --git a/boundary_conditions/dirichlets/dc_imp.py b/boundary_conditions/dirichlets/dc_imp.py deleted file mode 100644 index f53b43641..000000000 --- a/boundary_conditions/dirichlets/dc_imp.py +++ /dev/null @@ -1,95 +0,0 @@ -from festim import DirichletBC, BoundaryConditionExpression, k_B -import fenics as f -import sympy as sp - - -def dc_imp(T, phi, R_p, D_0, E_D, Kr_0=None, E_Kr=None, Kd_0=None, E_Kd=None, P=None): - D = D_0 * f.exp(-E_D / k_B / T) - value = phi * R_p / D - if Kr_0 is not None: - Kr = Kr_0 * f.exp(-E_Kr / k_B / T) - if Kd_0 is not None: - Kd = Kd_0 * f.exp(-E_Kd / k_B / T) - value += ((phi + Kd * P) / Kr) ** 0.5 - else: - value += (phi / Kr) ** 0.5 - - return value - - -class ImplantationDirichlet(DirichletBC): - """Subclass of DirichletBC representing an approximation of an implanted - flux of hydrogen. - The details of the approximation can be found in - https://www.nature.com/articles/s41598-020-74844-w - - c = phi*R_p/D + ((phi+Kd*P)/Kr)**0.5 - - Args: - surfaces (list or int): the surfaces of the BC - phi (float or sp.Expr): implanted flux (H/m2/s) - R_p (float or sp.Expr): implantation depth (m) - D_0 (float): diffusion coefficient pre-exponential factor (m2/s) - E_D (float): diffusion coefficient activation energy (eV) - Kr_0 (float, optional): recombination coefficient pre-exponential - factor (m^4/s). If None, instantaneous recombination will be - assumed. Defaults to None. - E_Kr (float, optional): recombination coefficient activation - energy (eV). Defaults to None. - Kd_0 (float, optional): dissociation coefficient pre-exponential - factor (m-2 s-1 Pa-1). If None, instantaneous dissociation will be - assumed. Defaults to None. - E_Kd (float, optional): dissociation coefficient activation - energy (eV). Defaults to None. - P (float or sp.Expr, optional): partial pressure of H (Pa). Defaults to None. - """ - - def __init__( - self, - surfaces, - phi, - R_p, - D_0, - E_D, - Kr_0=None, - E_Kr=None, - Kd_0=None, - E_Kd=None, - P=None, - ) -> None: - super().__init__(surfaces, field=0, value=None) - self.phi = phi - self.R_p = R_p - self.D_0 = D_0 - self.E_D = E_D - self.Kr_0 = Kr_0 - self.E_Kr = E_Kr - self.Kd_0 = Kd_0 - self.E_Kd = E_Kd - self.P = P - - def create_expression(self, T): - phi = f.Expression(sp.printing.ccode(self.phi), t=0, degree=1) - R_p = f.Expression(sp.printing.ccode(self.R_p), t=0, degree=1) - sub_expressions = [phi, R_p] - if self.P is not None: - P = f.Expression(sp.printing.ccode(self.P), t=0, degree=1) - sub_expressions.append(P) - else: - P = self.P - - value_BC = BoundaryConditionExpression( - T, - dc_imp, - phi=phi, - R_p=R_p, - D_0=self.D_0, - E_D=self.E_D, - Kr_0=self.Kr_0, - E_Kr=self.E_Kr, - Kd_0=self.Kd_0, - E_Kd=self.E_Kd, - P=P, - ) - self.expression = value_BC - self.sub_expressions = sub_expressions diff --git a/boundary_conditions/dirichlets/sieverts_bc.py b/boundary_conditions/dirichlets/sieverts_bc.py deleted file mode 100644 index 8920696b2..000000000 --- a/boundary_conditions/dirichlets/sieverts_bc.py +++ /dev/null @@ -1,37 +0,0 @@ -from festim import DirichletBC, BoundaryConditionExpression, k_B -import fenics as f -import sympy as sp - - -def sieverts_law(T, S_0, E_S, pressure): - S = S_0 * f.exp(-E_S / k_B / T) - return S * pressure**0.5 - - -class SievertsBC(DirichletBC): - """Subclass of DirichletBC for Sievert's law - - Args: - surfaces (list or int): the surfaces of the BC - S_0 (float): Sievert's constant pre-exponential factor (m-3/Pa0.5) - E_S (float): Sievert's constant activation energy (eV) - pressure (float or sp.Expr): hydrogen partial pressure (Pa) - """ - - def __init__(self, surfaces, S_0, E_S, pressure) -> None: - super().__init__(surfaces, field=0, value=None) - self.S_0 = S_0 - self.E_S = E_S - self.pressure = pressure - - def create_expression(self, T): - pressure = f.Expression(sp.printing.ccode(self.pressure), t=0, degree=1) - value_BC = BoundaryConditionExpression( - T, - sieverts_law, - S_0=self.S_0, - E_S=self.E_S, - pressure=pressure, - ) - self.expression = value_BC - self.sub_expressions = [pressure] diff --git a/boundary_conditions/fluxes/convective_flux.py b/boundary_conditions/fluxes/convective_flux.py deleted file mode 100644 index f81e49d98..000000000 --- a/boundary_conditions/fluxes/convective_flux.py +++ /dev/null @@ -1,27 +0,0 @@ -from festim import FluxBC -import fenics as f -import sympy as sp - - -class ConvectiveFlux(FluxBC): - """ - FluxBC subclass for convective heat flux - -lambda * grad(T) * n = h_coeff * (T - T_ext) - - Args: - h_coeff (float or sp.Expr): heat exchange coefficient (W/m2/K) - T_ext (float or sp.Expr): fluid temperature (K) - surfaces (list or int): the surfaces of the BC - """ - - def __init__(self, h_coeff, T_ext, surfaces) -> None: - self.h_coeff = h_coeff - self.T_ext = T_ext - super().__init__(surfaces=surfaces, field="T") - - def create_form(self, T, solute): - h_coeff = f.Expression(sp.printing.ccode(self.h_coeff), t=0, degree=1) - T_ext = f.Expression(sp.printing.ccode(self.T_ext), t=0, degree=1) - - self.form = -h_coeff * (T - T_ext) - self.sub_expressions = [h_coeff, T_ext] diff --git a/boundary_conditions/fluxes/dissociation_flux.py b/boundary_conditions/fluxes/dissociation_flux.py deleted file mode 100644 index 97954477a..000000000 --- a/boundary_conditions/fluxes/dissociation_flux.py +++ /dev/null @@ -1,33 +0,0 @@ -from festim import FluxBC, k_B -import fenics as f -import sympy as sp - - -class DissociationFlux(FluxBC): - """ - FluxBC subclass for hydrogen dissociation flux. - -D(T) * grad(c) * n = Kd(T) * P - - Args: - Kd_0 (float or sp.Expr): dissociation coefficient pre-exponential - factor (m-2 s-1 Pa-1) - E_Kd (float or sp.Expr): dissociation coefficient activation - energy (eV) - P (float or sp.Expr): partial pressure of H (Pa) - surfaces (list or int): the surfaces of the BC - """ - - def __init__(self, Kd_0, E_Kd, P, surfaces) -> None: - self.Kd_0 = Kd_0 - self.E_Kd = E_Kd - self.P = P - super().__init__(surfaces=surfaces, field=0) - - def create_form(self, T, solute): - Kd_0_expr = f.Expression(sp.printing.ccode(self.Kd_0), t=0, degree=1) - E_Kd_expr = f.Expression(sp.printing.ccode(self.E_Kd), t=0, degree=1) - P_expr = f.Expression(sp.printing.ccode(self.P), t=0, degree=1) - - Kd = Kd_0_expr * f.exp(-E_Kd_expr / k_B / T) - self.form = Kd * P_expr - self.sub_expressions = [Kd_0_expr, E_Kd_expr, P_expr] diff --git a/boundary_conditions/fluxes/mass_flux.py b/boundary_conditions/fluxes/mass_flux.py deleted file mode 100644 index 9e2e18fc5..000000000 --- a/boundary_conditions/fluxes/mass_flux.py +++ /dev/null @@ -1,30 +0,0 @@ -from festim import FluxBC -import fenics as f -import sympy as sp - - -class MassFlux(FluxBC): - """ - FluxBC subclass for advective mass flux - -D * grad(c) * n = h_mass * (c - c_ext) - - Args: - h_mass (float or sp.Expr): mass transfer coefficient (m/s) - c_ext (float or sp.Expr): external concentration (1/m3) - surfaces (list or int): the surfaces of the BC - - Reference: Bergman, T. L., Bergman, T. L., Incropera, F. P., Dewitt, D. P., - & Lavine, A. S. (2011). Fundamentals of heat and mass transfer. John Wiley & Sons. - """ - - def __init__(self, h_coeff, c_ext, surfaces) -> None: - self.h_coeff = h_coeff - self.c_ext = c_ext - super().__init__(surfaces=surfaces, field=0) - - def create_form(self, T, solute): - h_coeff = f.Expression(sp.printing.ccode(self.h_coeff), t=0, degree=1) - c_ext = f.Expression(sp.printing.ccode(self.c_ext), t=0, degree=1) - - self.form = -h_coeff * (solute - c_ext) - self.sub_expressions = [h_coeff, c_ext] diff --git a/boundary_conditions/fluxes/surface_kinetics.py b/boundary_conditions/fluxes/surface_kinetics.py deleted file mode 100644 index 53bf3f54f..000000000 --- a/boundary_conditions/fluxes/surface_kinetics.py +++ /dev/null @@ -1,184 +0,0 @@ -from festim import FluxBC -from fenics import * -import sympy as sp - - -class SurfaceKinetics(FluxBC): - r""" - FluxBC subclass allowing to include surface processes in 1D H transport simulations: - - .. math:: - - \dfrac{d c_{\mathrm{s}}}{dt} = J_{\mathrm{bs}} - J_{\mathrm{sb}} + J_{\mathrm{vs}}; - - .. math:: - - -D \nabla c_\mathrm{m} \cdot \mathbf{n} = \lambda_{\mathrm{IS}} \dfrac{\partial c_{\mathrm{m}}}{\partial t} - + J_{\mathrm{bs}} - J_{\mathrm{sb}}, - - where :math:`c_{\mathrm{m}}` is the concentration of mobile hydrogen (:math:`\mathrm{H \ m}^{-3}`), - :math:`c_{\mathrm{s}}` is the surface concentration of adsorbed hydrogen (:math:`\mathrm{H \ m}^{-2}`), - the H flux from subsurface to surface :math:`J_{\mathrm{bs}}` (in :math:`\mathrm{m}^{-2} \ \mathrm{s}^{-1}`) is: - - .. math:: - - J_{\mathrm{bs}} = k_{\mathrm{bs}} c_{\mathrm{m}} \lambda_{\mathrm{abs}} \left(1 - \dfrac{c_\mathrm{s}}{n_{\mathrm{surf}}}\right), - - the H flux from surface to subsurface :math:`J_{\mathrm{sb}}` (in :math:`\mathrm{m}^{-2} \ \mathrm{s}^{-1}`) is: - - .. math:: - - J_{\mathrm{sb}} = k_{\mathrm{sb}} c_{\mathrm{s}} \left(1 - \dfrac{c_{\mathrm{m}}}{n_\mathrm{IS}}\right), - - :math:`\lambda_{\mathrm{abs}}=n_{\mathrm{surf}}/n_{\mathrm{IS}}` is the characteristic distance between surface and - subsurface sites (:math:`\mathrm{m}`). - - For more details see: - E.A. Hodille et al 2017 Nucl. Fusion 57 056002; Y. Hamamoto et al 2020 Nucl. Mater. Energy 23 100751 - - .. warning:: - - The SurfaceKinetics boundary condition can be used only in 1D simulations! - - Args: - k_sb (float or callable): rate constant for the surface-to-subsurface transition (:math:`\mathrm{s}^{-1}`), - can accept additional parameters (see example) - k_bs (float or callable): rate constant for the subsurface-to-surface transition (:math:`\mathrm{s}^{-1}`), - can accept additional parameters (see example) - lambda_IS (float): characteristic distance between two iterstitial sites (:math:`\mathrm{m}`) - n_surf (float): surface concentration of adsorption sites (:math:`\mathrm{m}^{-2}`) - n_IS (float): bulk concentration of interstitial sites (:math:`\mathrm{m}^{-3}`) - J_vs (float or callable): the net adsorption flux from vacuum to surface (:math:`\mathrm{m}^{-2} \ \mathrm{s}^{-1}`), - can accept additional parameters (see example) - surfaces (int or list): the surfaces for which surface processes are considered - initial_condition (int or float): the initial value of the H surface concentration (:math:`\mathrm{m}^{-2}`) - - Attributes: - previous_solutions (list): list containing solutions (fenics.Function or ufl.Indexed) - on each surface for "previous" timestep - test_functions (list): list containing test functions (fenics.TestFunction or ufl.Indexed) - for each surface - post_processing_solutions (list): list containing solutions (fenics.Function or ufl.Indexed) - on each surface used for post-processing - - Example:: - - def K_sb(T, surf_conc, prm1, prm2): - return 1e13 * f.exp(-2.0/F.k_B/T) - - def K_bs(T, surf_conc, prm1, prm2): - return 1e13 * f.exp(-0.2/F.k_B/T) - - def J_vs(T, surf_conc, prm1, prm2): - return (1-surf_conc/5) ** 2 * fenics.exp(-2/F.k_B/T) + prm1 * prm2 - - my_surf_model = SurfaceKinetics( - k_sb=K_sb, - k_bs=K_bs, - lambda_IS=110e-12, - n_surf=2e19, - n_IS=6e28, - J_vs=J_vs, - surfaces=[1, 2], - initial_condition=0, - prm1=2e16, - prm2=F.t - ) - """ - - def __init__( - self, - k_sb, - k_bs, - lambda_IS, - n_surf, - n_IS, - J_vs, - surfaces, - initial_condition, - **prms, - ) -> None: - super().__init__(surfaces=surfaces, field=0) - self.k_sb = k_sb - self.k_bs = k_bs - self.J_vs = J_vs - self.lambda_IS = lambda_IS - self.n_surf = n_surf - self.n_IS = n_IS - self.J_vs = J_vs - self.initial_condition = initial_condition - self.prms = prms - self.convert_prms() - - self.solutions = [None] * len(self.surfaces) - self.previous_solutions = [None] * len(self.surfaces) - self.test_functions = [None] * len(self.surfaces) - self.post_processing_solutions = [None] * len(self.surfaces) - - def create_form(self, solute, solute_prev, solute_test_function, T, ds, dt): - """ - Creates the general form associated with the surface species - - Args: - solute (fenics.Function or ufl.Indexed): mobile solution for "current" - timestep - solute_prev (fenics.Function or ufl.Indexed): mobile solution for - "previous" timestep - solute_test_function (fenics.TestFunction or ufl.Indexed): mobile test function - T (festim.Temperature): the temperature of the simulation - ds (fenics.Measure): the ds measure of the sim - dt (festim.Stepsize): the step-size - """ - - lambda_IS = self.lambda_IS - n_surf = self.n_surf - n_IS = self.n_IS - lambda_abs = ( - n_surf / n_IS - ) # characteristic distance between surface and subsurface sites - self.form = 0 - - for i, surf in enumerate(self.surfaces): - - J_vs = self.J_vs - if callable(J_vs): - J_vs = J_vs(T.T, self.solutions[i], **self.prms) - k_sb = self.k_sb - if callable(k_sb): - k_sb = k_sb(T.T, self.solutions[i], **self.prms) - k_bs = self.k_bs - if callable(k_bs): - k_bs = k_bs(T.T, self.solutions[i], **self.prms) - - J_sb = k_sb * self.solutions[i] * (1 - solute / n_IS) - J_bs = k_bs * (solute * lambda_abs) * (1 - self.solutions[i] / n_surf) - - if dt is not None: - # Surface concentration form - self.form += ( - (self.solutions[i] - self.previous_solutions[i]) - / dt.value - * self.test_functions[i] - * ds(surf) - ) - # Flux to solute species - self.form += ( - lambda_IS - * (solute - solute_prev) - / dt.value - * solute_test_function - * ds(surf) - ) - - self.form += -(J_vs + J_bs - J_sb) * self.test_functions[i] * ds(surf) - self.form += (J_bs - J_sb) * solute_test_function * ds(surf) - - self.sub_expressions += [expression for expression in self.prms.values()] - - def convert_prms(self): - # create Expressions or Constant for all parameters - for key, value in self.prms.items(): - if isinstance(value, (int, float)): - self.prms[key] = Constant(value) - else: - self.prms[key] = Expression(sp.printing.ccode(value), t=0, degree=1) diff --git a/concentration/mobile.py b/concentration/mobile.py deleted file mode 100644 index ec2f666f1..000000000 --- a/concentration/mobile.py +++ /dev/null @@ -1,237 +0,0 @@ -from festim import Concentration, FluxBC, k_B, RadioactiveDecay, SurfaceKinetics -from fenics import * - - -class Mobile(Concentration): - """ - The mobile concentration. - - If conservation of chemical potential, this will be c_m/S. - If not, Mobile represents c_m. - - Attributes: - sources (list): list of festim.Source objects. - The volumetric source terms - F (fenics.Form): the variational formulation for mobile - """ - - def __init__(self): - """Inits festim.Mobile""" - super().__init__() - self.sources = [] - self.boundary_conditions = [] - - def create_form(self, materials, mesh, T, dt=None, traps=None, soret=False): - """Creates the variational formulation. - - Args: - materials (festim.Materials): the materials - mesh (festim.Mesh): the mesh of the simulation - T (festim.Temperature): the temperature - dt (festim.Stepsize, optional): the stepsize. Defaults to None. - traps (festim.Traps, optional): the traps. Defaults to None. - chemical_pot (bool, optional): if True, conservation of chemical - potential is assumed. Defaults to False. - soret (bool, optional): If True, Soret effect is assumed. Defaults - to False. - """ - self.F = 0 - self.create_diffusion_form(materials, mesh, T, dt=dt, traps=traps, soret=soret) - self.create_source_form(mesh.dx) - self.create_fluxes_form(T, mesh.ds, dt) - - def create_diffusion_form( - self, materials, mesh, T, dt=None, traps=None, soret=False - ): - """Creates the variational formulation for the diffusive part. - - Args: - materials (festim.Materials): the materials - mesh (festim.Mesh): the mesh - T (festim.Temperature): the temperature - dt (festim.Stepsize, optional): the stepsize. Defaults to None. - traps (festim.Traps, optional): the traps. Defaults to None. - chemical_pot (bool, optional): if True, conservation of chemical - potential is assumed. Defaults to False. - soret (bool, optional): If True, Soret effect is assumed. Defaults - to False. - """ - - F = 0 - for material in materials: - D_0 = material.D_0 - E_D = material.E_D - c_0, c_0_n = self.get_concentration_for_a_given_material(material, T) - - subdomains = material.id # list of subdomains with this material - if not isinstance(subdomains, list): - subdomains = [subdomains] # make sure subdomains is a list - - # add to the formulation F for every subdomain - for subdomain in subdomains: - dx = mesh.dx(subdomain) - # transient form - if dt is not None: - F += ((c_0 - c_0_n) / dt.value) * self.test_function * dx - D = D_0 * exp(-E_D / k_B / T.T) - if mesh.type == "cartesian": - F += dot(D * grad(c_0), grad(self.test_function)) * dx - if soret: - Q = material.Q - if callable(Q): - Q = Q(T.T) - F += ( - dot( - D * Q * c_0 / (k_B * T.T**2) * grad(T.T), - grad(self.test_function), - ) - * dx - ) - - # see https://fenicsproject.discourse.group/t/method-of-manufactured-solution-cylindrical/7963 - elif mesh.type == "cylindrical": - r = SpatialCoordinate(mesh.mesh)[0] - F += r * dot(D * grad(c_0), grad(self.test_function / r)) * dx - - if soret: - Q = material.Q - if callable(Q): - Q = Q(T.T) - F += ( - r - * dot( - D * Q * c_0 / (k_B * T.T**2) * grad(T.T), - grad(self.test_function / r), - ) - * dx - ) - - elif mesh.type == "spherical": - r = SpatialCoordinate(mesh.mesh)[0] - F += ( - D - * r - * r - * dot(grad(c_0), grad(self.test_function / r / r)) - * dx - ) - - if soret: - Q = material.Q - if callable(Q): - Q = Q(T.T) - F += ( - D - * r - * r - * dot( - Q * c_0 / (k_B * T.T**2) * grad(T.T), - grad(self.test_function / r / r), - ) - * dx - ) - - # add the trapping terms - F_trapping = 0 - if traps is not None: - for trap in traps: - for i, mat in enumerate(trap.materials): - if type(trap.k_0) is list: - k_0 = trap.k_0[i] - E_k = trap.E_k[i] - p_0 = trap.p_0[i] - E_p = trap.E_p[i] - density = trap.density[i] - else: - k_0 = trap.k_0 - E_k = trap.E_k - p_0 = trap.p_0 - E_p = trap.E_p - density = trap.density[0] - c_m, _ = self.get_concentration_for_a_given_material(mat, T) - F_trapping += ( - -k_0 - * exp(-E_k / k_B / T.T) - * c_m - * (density - trap.solution) - * self.test_function - * dx(mat.id) - ) - F_trapping += ( - p_0 - * exp(-E_p / k_B / T.T) - * trap.solution - * self.test_function - * dx(mat.id) - ) - F += -F_trapping - - self.F_diffusion = F - self.F += F - - def create_source_form(self, dx): - """Creates the variational form for the volumetric source term parts. - - Args: - dx (fenics.Measure): the measure dx - """ - F_source = 0 - expressions_source = [] - - print("Defining source terms") - for source in self.sources: - if type(source.volume) is list: - volumes = source.volume - else: - volumes = [source.volume] - if isinstance(source, RadioactiveDecay): - source.value = source.form(self.mobile_concentration()) - - for volume in volumes: - F_source += -source.value * self.test_function * dx(volume) - if isinstance(source.value, (Expression, UserExpression)): - expressions_source.append(source.value) - - self.F_source = F_source - self.F += F_source - self.sub_expressions += expressions_source - - def create_fluxes_form(self, T, ds, dt=None): - """Modifies the formulation and adds fluxes based - on parameters in self.boundary_conditions - """ - - expressions_fluxes = [] - F = 0 - - solute = self.mobile_concentration() - - for bc in self.boundary_conditions: - if bc.field != "T": - if isinstance(bc, FluxBC): - if isinstance(bc, SurfaceKinetics): - bc.create_form( - solute, - self.previous_solution, - self.test_function, - T, - ds, - dt, - ) - F += bc.form - else: - bc.create_form(T.T, solute) - for surf in bc.surfaces: - F += -self.test_function * bc.form * ds(surf) - # TODO : one day we will get rid of this huge expressions list - expressions_fluxes += bc.sub_expressions - - self.F_fluxes = F - self.F += F - self.sub_expressions += expressions_fluxes - - def get_concentration_for_a_given_material(self, material, T): - return self.solution, self.previous_solution - - def mobile_concentration(self): - return self.solution diff --git a/concentration/theta.py b/concentration/theta.py deleted file mode 100644 index 51fb81615..000000000 --- a/concentration/theta.py +++ /dev/null @@ -1,122 +0,0 @@ -from festim import Mobile, k_B -import fenics as f - - -class Theta(Mobile): - """Class representing the "chemical potential" c/S where S is the - solubility of the metal - """ - - def __init__(self): - """Inits Theta""" - super().__init__() - self.S = None - self.F = None - - def initialise(self, V, value, label=None, time_step=None): - """Assign a value to self.previous_solution - - Args: - V (fenics.FunctionSpace): the function space - value (sp.Add, float, int, str): the value of the initialisation. - label (str, optional): the label in the XDMF file. Defaults to - None. - time_step (int, optional): the time step to read in the XDMF file. - Defaults to None. - """ - comp = self.get_comp(V, value, label=label, time_step=time_step) - - prev_sol = f.Function(V) - v = f.TestFunction(V) - dx = f.Measure("dx", subdomain_data=self.volume_markers) - F = 0 - for mat in self.materials: - S = mat.S_0 * f.exp(-mat.E_S / k_B / self.T.T) - F += -prev_sol * v * dx(mat.id) - if mat.solubility_law == "sievert": - F += comp / S * v * dx(mat.id) - elif mat.solubility_law == "henry": - F += (comp / S) ** 0.5 * v * dx(mat.id) - f.solve(F == 0, prev_sol, bcs=[]) - - f.assign(self.previous_solution, prev_sol) - - def get_concentration_for_a_given_material(self, material, T): - """Returns the concentration (and previous concentration) for a given - material - - Args: - material (festim.Material): the material with attributes S_0 and - E_S - T (festim.Temperature): the temperature with attributest T and T_n - - Returns: - fenics.Product, fenics.Product: the current concentration and - previous concentration - """ - E_S = material.E_S - S_0 = material.S_0 - S = S_0 * f.exp(-E_S / k_B / T.T) - S_n = S_0 * f.exp(-E_S / k_B / T.T_n) - if material.solubility_law == "sievert": - c_0 = self.solution * S - c_0_n = self.previous_solution * S_n - elif material.solubility_law == "henry": - c_0 = (self.solution) ** 2 * S - c_0_n = self.previous_solution**2 * S_n - return c_0, c_0_n - - def mobile_concentration(self): - """Returns the hydrogen concentration as c=theta*K_S or c=theta**2*K_H - This is needed when adding robin BCs (eg RecombinationFlux). - - Returns: - ufl.algebra.Sum: the hydrogen mobile concentration - """ - henry_to_concentration = self.solution**2 * self.S - sieverts_to_concentration = self.solution * self.S - # henry_marker is equal to 1 in Henry materials and 0 elsewhere - return ( - self.materials.henry_marker * henry_to_concentration - + self.materials.sievert_marker * sieverts_to_concentration - ) - - def post_processing_solution_to_concentration(self): - """Converts the post_processing_solution from theta to mobile - concentration. - c = theta * S. - The attribute post_processing_solution is fenics.Product (if self.S is - festim.ArheniusCoeff) - """ - problem = f.LinearVariationalProblem( - a=f.lhs(self.form_post_processing), - L=f.rhs(self.form_post_processing), - u=self.post_processing_solution, - bcs=[], - ) - solver = f.LinearVariationalSolver(problem) - solver.solve() - - def create_form_post_processing(self, V, materials, dx): - """Creates a variational formulation for c = theta * S or theta**2 * S - - Args: - V (fenics.FunctionSpace): the DG1 function space of the concentration field - materials (festim.Materials): the materials - dx (fenics.Measurement): the dx measure of the problem - """ - F = 0 - v = f.TestFunction(V) - c = f.TrialFunction(V) - - F += -c * v * dx - for mat in materials: - if mat.solubility_law == "sievert": - # for sievert materials c = theta * S - F += self.solution * self.S * v * dx(mat.id) - elif mat.solubility_law == "henry": - # for henry materials c = theta**2 * S - F += self.solution**2 * self.S * v * dx(mat.id) - - self.form_post_processing = F - self.post_processing_solution = f.Function(V) diff --git a/concentration/traps/extrinsic_trap.py b/concentration/traps/extrinsic_trap.py deleted file mode 100644 index 78a450309..000000000 --- a/concentration/traps/extrinsic_trap.py +++ /dev/null @@ -1,161 +0,0 @@ -from festim import Trap, as_constant_or_expression -from fenics import NewtonSolver, MPI - - -class ExtrinsicTrapBase(Trap): - def __init__( - self, - k_0, - E_k, - p_0, - E_p, - materials, - id=None, - absolute_tolerance=1e0, - relative_tolerance=1e-10, - maximum_iterations=30, - linear_solver=None, - preconditioner="default", - **kwargs, - ): - """Inits ExtrinsicTrap - - Args: - E_k (float, list): trapping pre-exponential factor (m3 s-1) - k_0 (float, list): trapping activation energy (eV) - p_0 (float, list): detrapping pre-exponential factor (s-1) - E_p (float, list): detrapping activation energy (eV) - materials (list, int): the materials ids the trap is living in - id (int, optional): The trap id. Defaults to None. - absolute_tolerance (float, optional): the absolute tolerance of the newton - solver. Defaults to 1e-0 - relative_tolerance (float, optional): the relative tolerance of the newton - solver. Defaults to 1e-10 - maximum_iterations (int, optional): maximum iterations allowed for - the solver to converge. Defaults to 30. - linear_solver (str, optional): linear solver method for the newton solver, - options can be viewed with print(list_linear_solver_methods()). - If None, the default fenics linear solver will be used ("umfpack"). - More information can be found at: https://fenicsproject.org/pub/tutorial/html/._ftut1017.html. - Defaults to None. - preconditioner (str, optional): preconditioning method for the newton solver, - options can be viewed by print(list_krylov_solver_preconditioners()). - Defaults to "default". - """ - super().__init__(k_0, E_k, p_0, E_p, materials, density=None, id=id) - self.absolute_tolerance = absolute_tolerance - self.relative_tolerance = relative_tolerance - self.maximum_iterations = maximum_iterations - self.linear_solver = linear_solver - self.preconditioner = preconditioner - - self.newton_solver = None - for name, val in kwargs.items(): - setattr(self, name, as_constant_or_expression(val)) - self.density_previous_solution = None - self.density_test_function = None - - @property - def newton_solver(self): - return self._newton_solver - - @newton_solver.setter - def newton_solver(self, value): - if value is None: - self._newton_solver = value - elif isinstance(value, NewtonSolver): - if self._newton_solver: - print("Settings for the Newton solver will be overwritten") - self._newton_solver = value - else: - raise TypeError("accepted type for newton_solver is fenics.NewtonSolver") - - def define_newton_solver(self): - """Creates the Newton solver and sets its parameters""" - self.newton_solver = NewtonSolver(MPI.comm_world) - self.newton_solver.parameters["error_on_nonconvergence"] = True - self.newton_solver.parameters["absolute_tolerance"] = self.absolute_tolerance - self.newton_solver.parameters["relative_tolerance"] = self.relative_tolerance - self.newton_solver.parameters["maximum_iterations"] = self.maximum_iterations - self.newton_solver.parameters["linear_solver"] = self.linear_solver - self.newton_solver.parameters["preconditioner"] = self.preconditioner - - -class ExtrinsicTrap(ExtrinsicTrapBase): - """ - For details in the forumation see - http://www.sciencedirect.com/science/article/pii/S2352179119300547 - - Args: - E_k (float, list): trapping pre-exponential factor (m3 s-1) - k_0 (float, list): trapping activation energy (eV) - p_0 (float, list): detrapping pre-exponential factor (s-1) - E_p (float, list): detrapping activation energy (eV) - materials (list, int): the materials ids the trap is living in - id (int, optional): The trap id. Defaults to None. - """ - - def __init__( - self, - k_0, - E_k, - p_0, - E_p, - materials, - phi_0, - n_amax, - n_bmax, - eta_a, - eta_b, - f_a, - f_b, - id=None, - **kwargs, - ): - super().__init__( - k_0, - E_k, - p_0, - E_p, - materials, - phi_0=phi_0, - n_amax=n_amax, - n_bmax=n_bmax, - eta_a=eta_a, - eta_b=eta_b, - f_a=f_a, - f_b=f_b, - id=id, - **kwargs, - ) - - def create_form_density(self, dx, dt, T): - """ - Creates the variational formulation for the extrinsic trap density. - - Args: - dx (fenics.Measure): the dx measure of the sim - dt (festim.Stepsize): the stepsize of the simulation. - T (festim.Temperature): the temperature of the - simulation - - .. note:: - T is an argument, although is not used in the formulation of - extrinsic traps, but potential for subclasses of extrinsic traps - """ - density = self.density[0] - F = ( - ((density - self.density_previous_solution) / dt.value) - * self.density_test_function - * dx - ) - F += ( - -self.phi_0 - * ( - (1 - density / self.n_amax) * self.eta_a * self.f_a - + (1 - density / self.n_bmax) * self.eta_b * self.f_b - ) - * self.density_test_function - * dx - ) - self.form_density = F diff --git a/concentration/traps/trap.py b/concentration/traps/trap.py deleted file mode 100644 index 1fef133ba..000000000 --- a/concentration/traps/trap.py +++ /dev/null @@ -1,211 +0,0 @@ -from festim import Concentration, k_B, Material, Theta, RadioactiveDecay -from fenics import * -import sympy as sp - - -class Trap(Concentration): - """ - Args: - k_0 (float, list): trapping pre-exponential factor (m3 s-1) - E_k (float, list): trapping activation energy (eV) - p_0 (float, list): detrapping pre-exponential factor (s-1) - E_p (float, list): detrapping activation energy (eV) - materials (list, str, festim.Material): the materials the - trap is living in. The material's name. - density (sp.Add, float, list, fenics.Expresion, fenics.UserExpression): - the trap density (m-3) - id (int, optional): The trap id. Defaults to None. - - Raises: - ValueError: if duplicates are found in materials - - .. note:: - Should multiple traps in muliple materials be used, to save on - dof's, traps can be conglomerated and described in lists in the - format:: - - festim.Trap( - k_0=[1, 2], - E_k=[1, 2], - p_0=[1, 2], - E_p=[1, 2], - materials=[1, 2] - density=[1, 2]) - - This will act as a singular trap but with seperate properties for - respective materials. Parameters k_0, E_k, p_0, E_p, materials and - density MUST have the same length for this method to be valid. - """ - - def __init__(self, k_0, E_k, p_0, E_p, materials, density, id=None): - super().__init__() - self.id = id - self.k_0 = k_0 - self.E_k = E_k - self.p_0 = p_0 - self.E_p = E_p - self.materials = materials - - self.density = [] - self.make_density(density) - self.sources = [] - - @property - def materials(self): - return self._materials - - @materials.setter - def materials(self, value): - if not isinstance(value, list): - value = [value] - for entry in value: - if not isinstance(entry, (str, Material)): - raise TypeError( - "Accepted types for materials are str or festim.Material" - ) - self._materials = value - - def make_materials(self, materials): - """Ensure all entries in self.materials are of type festim.Material - - Args: - materials (festim.Materials): the materials - - Raises: - ValueError: if some duplicates are found - """ - new_materials = [] - - for material in self.materials: - new_materials.append(materials.find_material(material)) - - self.materials = new_materials - - if len(self.materials) != len(list(set(self.materials))): - raise ValueError("Duplicate materials in trap") - - def make_density(self, densities): - if type(densities) is not list: - densities = [densities] - - for i, density in enumerate(densities): - if density is not None: - # if density is already a fenics Expression, use it as is - if isinstance(density, (Expression, UserExpression)): - self.density.append(density) - # else assume it's a sympy expression - else: - density_expr = sp.printing.ccode(density) - self.density.append( - Expression( - density_expr, - degree=2, - t=0, - name="density_{}_{}".format(self.id, i), - ) - ) - - def create_form(self, mobile, materials, T, dx, dt=None): - """Creates the general form associated with the trap - d ct/ dt = k c_m (n - c_t) - p c_t + S - - Args: - mobile (festim.Mobile): the mobile concentration of the simulation - materials (festim.Materials): the materials of the simulation - T (festim.Temperature): the temperature of the simulation - dx (fenics.Measure): the dx measure of the sim - dt (festim.Stepsize, optional): If None assuming steady state. - Defaults to None. - """ - self.F = 0 - self.create_trapping_form(mobile, materials, T, dx, dt) - if self.sources is not None: - self.create_source_form(dx) - - def create_trapping_form(self, mobile, materials, T, dx, dt=None): - """d ct/ dt = k c_m (n - c_t) - p c_t - - Args: - mobile (festim.Mobile): the mobile concentration of the simulation - materials (festim.Materials): the materials of the simulation - T (festim.Temperature): the temperature of the simulation - dx (fenics.Measure): the dx measure of the sim - dt (festim.Stepsize, optional): If None assuming steady state. - Defaults to None. - """ - solution = self.solution - prev_solution = self.previous_solution - test_function = self.test_function - - if not all(isinstance(mat, Material) for mat in self.materials): - self.make_materials(materials) - - expressions_trap = [] - F_trapping = 0 # initialise the form - - if dt is not None: - # d(c_t)/dt in trapping equation - F_trapping += ((solution - prev_solution) / dt.value) * test_function * dx - else: - # if the sim is steady state and - # if a trap is not defined in one subdomain - # add c_t = 0 to the form in this subdomain - for mat in materials: - if mat not in self.materials: - F_trapping += solution * test_function * dx(mat.id) - - for i, mat in enumerate(self.materials): - if type(self.k_0) is list: - k_0 = self.k_0[i] - E_k = self.E_k[i] - p_0 = self.p_0[i] - E_p = self.E_p[i] - density = self.density[i] - else: - k_0 = self.k_0 - E_k = self.E_k - p_0 = self.p_0 - E_p = self.E_p - density = self.density[0] - - # add the density to the list of - # expressions to be updated - expressions_trap.append(density) - - if isinstance(mobile, Theta) and mat.solubility_law == "henry": - raise NotImplementedError( - "Henry law of solubility is not implemented with traps" - ) - - c_0, c_0_n = mobile.get_concentration_for_a_given_material(mat, T) - - # k(T)*c_m*(n - c_t) - p(T)*c_t - F_trapping += ( - -k_0 - * exp(-E_k / k_B / T.T) - * c_0 - * (density - solution) - * test_function - * dx(mat.id) - ) - F_trapping += ( - p_0 * exp(-E_p / k_B / T.T) * solution * test_function * dx(mat.id) - ) - - self.F_trapping = F_trapping - self.F += self.F_trapping - self.sub_expressions += expressions_trap - - def create_source_form(self, dx): - """Create the source form for the trap - - Args: - dx (fenics.Measure): the dx measure of the sim - """ - for source in self.sources: - if isinstance(source, RadioactiveDecay): - source.value = source.form(self.solution) - self.F_source = -source.value * self.test_function * dx(source.volume) - self.F += self.F_source - if isinstance(source.value, (Expression, UserExpression)): - self.sub_expressions.append(source.value) diff --git a/concentration/traps/traps.py b/concentration/traps/traps.py deleted file mode 100644 index e819ff4cf..000000000 --- a/concentration/traps/traps.py +++ /dev/null @@ -1,133 +0,0 @@ -import festim -import fenics as f -import warnings - - -class Traps(list): - """ - A list of festim.Trap objects - """ - - def __init__(self, *args): - # checks that input is list - if len(args) == 0: - super().__init__() - else: - if not isinstance(*args, list): - raise TypeError("festim.Traps must be a list") - super().__init__(self._validate_trap(item) for item in args[0]) - - self.F = None - self.extrinsic_formulations = [] - self.sub_expressions = [] - - @property - def traps(self): - warnings.warn( - "The traps attribute will be deprecated in a future release, please use festim.Traps as a list instead", - DeprecationWarning, - ) - return self - - @traps.setter - def traps(self, value): - warnings.warn( - "The traps attribute will be deprecated in a future release, please use festim.Traps as a list instead", - DeprecationWarning, - ) - if isinstance(value, list): - if not all(isinstance(t, festim.Trap) for t in value): - raise TypeError("traps must be a list of festim.Trap") - super().__init__(value) - else: - raise TypeError("traps must be a list") - - def __setitem__(self, index, item): - super().__setitem__(index, self._validate_trap(item)) - - def insert(self, index, item): - super().insert(index, self._validate_trap(item)) - - def append(self, item): - super().append(self._validate_trap(item)) - - def extend(self, other): - if isinstance(other, type(self)): - super().extend(other) - else: - super().extend(self._validate_trap(item) for item in other) - - def _validate_trap(self, value): - if isinstance(value, festim.Trap): - return value - raise TypeError("festim.Traps must be a list of festim.Trap") - - def make_traps_materials(self, materials): - for trap in self: - trap.make_materials(materials) - - def assign_traps_ids(self): - for i, trap in enumerate(self, 1): - if trap.id is None: - trap.id = i - - def create_forms(self, mobile, materials, T, dx, dt=None): - self.F = 0 - for trap in self: - trap.create_form(mobile, materials, T, dx, dt=dt) - self.F += trap.F - self.sub_expressions += trap.sub_expressions - - def get_trap(self, id): - for trap in self: - if trap.id == id: - return trap - raise ValueError("Couldn't find trap {}".format(id)) - - def initialise_extrinsic_traps(self, V): - """Add functions to ExtrinsicTrapBase objects for density form""" - for trap in self: - if isinstance(trap, festim.ExtrinsicTrapBase): - trap.density = [f.Function(V)] - trap.density_test_function = f.TestFunction(V) - trap.density_previous_solution = f.project(f.Constant(0), V) - - def define_variational_problem_extrinsic_traps(self, dx, dt, T): - """ - Creates the variational formulations for the extrinsic traps densities - - Args: - dx (fenics.Measure): the dx measure of the sim - dt (festim.Stepsize): If None assuming steady state. - T (festim.Temperature): the temperature of the simulation - """ - self.extrinsic_formulations = [] - expressions_extrinsic = [] - for trap in self: - if isinstance(trap, festim.ExtrinsicTrapBase): - trap.create_form_density(dx, dt, T) - self.extrinsic_formulations.append(trap.form_density) - self.sub_expressions.extend(expressions_extrinsic) - - def define_newton_solver_extrinsic_traps(self): - for trap in self: - if isinstance(trap, festim.ExtrinsicTrapBase): - trap.define_newton_solver() - - def solve_extrinsic_traps(self): - for trap in self: - if isinstance(trap, festim.ExtrinsicTrapBase): - du_t = f.TrialFunction(trap.density[0].function_space()) - J_t = f.derivative(trap.form_density, trap.density[0], du_t) - problem = festim.Problem(J_t, trap.form_density, []) - - f.begin( - "Solving nonlinear variational problem." - ) # Add message to fenics logs - trap.newton_solver.solve(problem, trap.density[0].vector()) - f.end() - - def update_extrinsic_traps_density(self): - for trap in self: - if isinstance(trap, festim.ExtrinsicTrapBase): - trap.density_previous_solution.assign(trap.density[0]) diff --git a/exports/derived_quantities/adsorbed_hydrogen.py b/exports/derived_quantities/adsorbed_hydrogen.py deleted file mode 100644 index 4f9f5d076..000000000 --- a/exports/derived_quantities/adsorbed_hydrogen.py +++ /dev/null @@ -1,43 +0,0 @@ -from festim import SurfaceQuantity -import fenics as f - - -class AdsorbedHydrogen(SurfaceQuantity): - """ - Object to compute the value of the adsorbed H concentration, - defined with the ``SurfaceKinetics`` boundary condition on a given surface. - - .. warning:: - - The ``AdsorbedHydrogen`` export can be used only if the ``SurfaceKinetics`` - condition is defined on the same surface! - - Args: - surface (int): the surface id - - Attributes: - surface (int): the surface id - export_unit (str): the unit of the derived quantity in the export file - title (str): the title of the derived quantity - show_units (bool): show the units in the title in the derived quantities file - function (dolfin.function.function.Function): the solution function of the hydrogen adsorbed field - - """ - - def __init__(self, surface) -> None: - super().__init__(field="adsorbed", surface=surface) - - @property - def export_unit(self): - return f"H m-2" - - @property - def title(self): - quantity_title = f"Adsorbed H on surface {self.surface}" - if self.show_units: - return quantity_title + f" ({self.export_unit})" - else: - return quantity_title - - def compute(self): - return f.assemble(self.function * self.ds(self.surface)) diff --git a/exports/derived_quantities/average_surface.py b/exports/derived_quantities/average_surface.py deleted file mode 100644 index c45f67218..000000000 --- a/exports/derived_quantities/average_surface.py +++ /dev/null @@ -1,49 +0,0 @@ -from festim import SurfaceQuantity -import fenics as f - - -class AverageSurface(SurfaceQuantity): - """ - Computes the average value of a field on a given surface - int(f ds) / int (1 * ds) - - Args: - field (str, int): the field ("solute", 0, 1, "T", "retention") - surface (int): the surface id - - Attributes: - field (str, int): the field ("solute", 0, 1, "T", "retention") - surface (int): the surface id - title (str): the title of the derived quantity - show_units (bool): show the units in the title in the derived quantities - file - function (dolfin.function.function.Function): the solution function of - the field - - .. note:: - Units are in H/m3 for hydrogen concentration and K for temperature - - """ - - def __init__(self, field, surface) -> None: - super().__init__(field=field, surface=surface) - - @property - def allowed_meshes(self): - return ["cartesian"] - - @property - def title(self): - quantity_title = f"Average {self.field} surface {self.surface}" - if self.show_units: - if self.field == "T": - return quantity_title + " (K)" - else: - return quantity_title + " (H m-3)" - else: - return quantity_title - - def compute(self): - return f.assemble(self.function * self.ds(self.surface)) / f.assemble( - 1 * self.ds(self.surface) - ) diff --git a/exports/derived_quantities/average_volume.py b/exports/derived_quantities/average_volume.py deleted file mode 100644 index 7bc0e07c2..000000000 --- a/exports/derived_quantities/average_volume.py +++ /dev/null @@ -1,48 +0,0 @@ -from festim import VolumeQuantity -import fenics as f - - -class AverageVolume(VolumeQuantity): - """ - Computes the average value of a field in a given volume - int(f dx) / int (1 * dx) - - Args: - field (str, int): the field ("solute", 0, 1, "T", "retention") - volume (int): the volume id - - Attributes: - field (str, int): the field ("solute", 0, 1, "T", "retention") - volume (int): the volume id - title (str): the title of the derived quantity - show_units (bool): show the units in the title in the derived quantities - file - function (dolfin.function.function.Function): the solution function of - the field - - .. note:: - Units are in H/m3 for hydrogen concentration and K for temperature - """ - - def __init__(self, field, volume: int) -> None: - super().__init__(field=field, volume=volume) - - @property - def allowed_meshes(self): - return ["cartesian"] - - @property - def title(self): - quantity_title = f"Average {self.field} volume {self.volume}" - if self.show_units: - if self.field == "T": - return quantity_title + " (K)" - else: - return quantity_title + " (H m-3)" - else: - return quantity_title - - def compute(self): - return f.assemble(self.function * self.dx(self.volume)) / f.assemble( - 1 * self.dx(self.volume) - ) diff --git a/exports/derived_quantities/derived_quantities.py b/exports/derived_quantities/derived_quantities.py deleted file mode 100644 index 28a8a3608..000000000 --- a/exports/derived_quantities/derived_quantities.py +++ /dev/null @@ -1,302 +0,0 @@ -from festim import ( - MinimumVolume, - MaximumVolume, - DerivedQuantity, -) -import fenics as f -import os -import numpy as np -from typing import Union -import warnings - - -class DerivedQuantities(list): - """ - A list of festim.DerivedQuantity objects - - Args: - filename (str, optional): the filename (must end with .csv). - If None, the data will not be exported. Defaults to None. - nb_iterations_between_compute (int, optional): number of - iterations between each derived quantities computation. - Defaults to 1. - nb_iterations_between_exports (int, optional): number of - iterations between each export. If None, the file will be - exported at the last timestep. Defaults to None. - show_units (bool, optional): will show the units of each - derived quantity in the title in export - - Attributes: - filename (str): the filename. - nb_iterations_between_compute (int): number of - iterations between each derived quantities computation. - nb_iterations_between_exports (int): number of - iterations between each export. If None, the file will be - exported at the last timestep. - show_units (bool): will show the units of each - derived quantity in the title in export - data (list): the data to be exported - t (list): the time steps - """ - - def __init__( - self, - *args, - filename: str = None, - nb_iterations_between_compute: int = 1, - nb_iterations_between_exports: int = None, - show_units=False, - ) -> None: - # checks that input is list - if len(args) == 0: - super().__init__() - else: - if not isinstance(*args, list): - raise TypeError("festim.DerivedQuantities must be a list") - super().__init__(self._validate_derived_quantity(item) for item in args[0]) - - self.filename = filename - self.nb_iterations_between_compute = nb_iterations_between_compute - self.nb_iterations_between_exports = nb_iterations_between_exports - self.show_units = show_units - - self.data = [] - self.t = [] - - @property - def derived_quantities(self): - warnings.warn( - "The derived_quantities attribute will be deprecated in a future release, please use festim.DerivedQuantities as a list instead", - DeprecationWarning, - ) - return self - - @derived_quantities.setter - def derived_quantities(self, value): - warnings.warn( - "The derived_quantities attribute will be deprecated in a future release, please use festim.DerivedQuantities as a list instead", - DeprecationWarning, - ) - if isinstance(value, list): - if not all(isinstance(t, DerivedQuantity) for t in value): - raise TypeError( - "derived_quantities must be a list of festim.DerivedQuantity" - ) - super().__init__(value) - else: - raise TypeError("derived_quantities must be a list") - - def __setitem__(self, index, item): - super().__setitem__(index, self._validate_derived_quantity(item)) - - def insert(self, index, item): - super().insert(index, self._validate_derived_quantity(item)) - - def append(self, item): - super().append(self._validate_derived_quantity(item)) - - def extend(self, other): - if isinstance(other, type(self)): - super().extend(other) - else: - super().extend(self._validate_derived_quantity(item) for item in other) - - def _validate_derived_quantity(self, value): - if isinstance(value, DerivedQuantity): - return value - raise TypeError( - "festim.DerivedQuantities must be a list of festim.DerivedQuantity" - ) - - @property - def filename(self): - return self._filename - - @filename.setter - def filename(self, value): - if value is not None: - if not isinstance(value, str): - raise TypeError("filename must be a string") - if not value.endswith(".csv"): - raise ValueError("filename must end with .csv") - self._filename = value - - def make_header(self): - header = ["t(s)"] - for quantity in self: - quantity.show_units = self.show_units - if self.show_units is False: - warnings.warn( - "The current derived_quantities title style will be deprecated in a future release, please use show_units=True instead", - DeprecationWarning, - ) - header.append(quantity.title) - return header - - def assign_measures_to_quantities(self, dx, ds): - self.volume_markers = dx.subdomain_data() - for quantity in self: - quantity.dx = dx - quantity.ds = ds - quantity.n = f.FacetNormal(dx.subdomain_data().mesh()) - - def assign_properties_to_quantities(self, materials): - """Assign properties attributes to all DerivedQuantity objects - (D, S, thermal_cond and H) based on the properties stored in materials - - Args: - materials (festim.Materials): the materials - """ - for quantity in self: - quantity.D = materials.D - quantity.S = materials.S - quantity.thermal_cond = materials.thermal_cond - quantity.Q = materials.Q - - def compute(self, t): - row = [t] - for quantity in self: - if isinstance(quantity, (MaximumVolume, MinimumVolume)): - value = quantity.compute(self.volume_markers) - else: - value = quantity.compute() - - # check if first time writing data - if len(self.data) == 0: - self.data = [self.make_header()] - - quantity.data.append(value) - quantity.t.append(t) - row.append(value) - self.data.append(row) - self.t.append(t) - - def write(self): - if self.filename is not None: - # 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) - - # save the data to csv - np.savetxt(self.filename, np.array(self.data), fmt="%s", delimiter=",") - return True - - def is_export(self, t, final_time, nb_iterations): - """Checks if the derived quantities should be exported or not based on - the current time, the final time of simulation and the current number - of iterations - - Args: - t (float): the current time - final_time (float): the final time of the simulation - nb_iterations (int): the current number of time steps - - Returns: - bool: True if the derived quantities should be exported, else False - """ - if final_time is not None: - nb_its_between_exports = self.nb_iterations_between_exports - if nb_its_between_exports is None: - # export at the end - return np.isclose(t, final_time, atol=0) - else: - # export every N iterations - return nb_iterations % nb_its_between_exports == 0 - else: - # if steady state, export - return True - - def is_compute(self, nb_iterations): - """Checks if the derived quantities should be computed or not based on - the current number of iterations - - Args: - nb_iterations (int): the current number of time steps - - Returns: - bool: True if it's time to compute, else False - """ - return nb_iterations % self.nb_iterations_between_compute == 0 - - def filter( - self, - surfaces: Union[list, int] = None, - volumes: Union[list, int] = None, - fields: Union[list, str] = None, - instances: DerivedQuantity = None, - ): - """Finds DerivedQuantity objects that match surfaces, volumes, and instances. - - Args: - surfaces (Union[list, int], optional): the surface ids to match. - Defaults to None. - volumes (Union[list, int], optional): the volume ids to match. - Defaults to None. - fields (Union[list, str], optional): the fields to match. - Defaults to None. - instances (DerivedQuantity, optional): the DerivedQuantity - instances to match. Defaults to None. - - Returns: - list, DerivedQuantity: if only one quantity matches returns this - quantity, else returs a list of DerivedQuantity - """ - # ensure arguments are list - if surfaces is not None and not isinstance(surfaces, list): - surfaces = [surfaces] - if volumes is not None and not isinstance(volumes, list): - volumes = [volumes] - if fields is not None and not isinstance(fields, list): - fields = [fields] - if instances is not None and not isinstance(instances, list): - instances = [instances] - - quantities = [] - - # iterate through derived_quantities - for quantity in self: - # initialise flags to False - match_surface, match_volume, match_field, match_instance = ( - False, - False, - False, - False, - ) - - # check if matches surface - if surfaces is not None: - if hasattr(quantity, "surface") and quantity.surface in surfaces: - match_surface = True - else: - match_surface = True - - # check if matches volume - if volumes is not None: - if hasattr(quantity, "volume") and quantity.volume in volumes: - match_volume = True - else: - match_volume = True - - # check if matches field - if fields is not None: - if quantity.field in fields: - match_field = True - else: - match_field = True - - # check if matches instance - if instances is not None: - if isinstance(quantity, tuple(instances)): - match_instance = True - else: - match_instance = True - - # if all flags are True, append to the list - if match_surface and match_volume and match_field and match_instance: - quantities.append(quantity) - - if len(quantities) == 1: - quantities = quantities[0] - return quantities diff --git a/exports/derived_quantities/derived_quantity.py b/exports/derived_quantities/derived_quantity.py deleted file mode 100644 index e4bf5a1b8..000000000 --- a/exports/derived_quantities/derived_quantity.py +++ /dev/null @@ -1,98 +0,0 @@ -from festim import Export - - -class DerivedQuantity(Export): - """ - Parent class of all derived quantities - - Args: - field (str, int): the field ("solute", 0, 1, "T", "retention") - - Attributes: - field (str, int): the field ("solute", 0, 1, "T", "retention") - title (str): the title of the derived quantity - show_units (bool): show the units in the title in the derived quantities - file - function (fenics.Function): the solution function of - the field - dx (fenics.Measure): the measure of the volume - ds (fenics.Measure): the measure of the surface - n (fenics.Function): the normal vector - D (fenics.Function): the diffusion coefficient - S (fenics.Function): the source term - thermal_cond (fenics.Function): the thermal conductivity - Q (fenics.Function): the heat source term - data (list): the data of the derived quantity - t (list): the time values of the data - allowed_meshes (list): the allowed meshes for the derived quantity - """ - - def __init__(self, field) -> None: - super().__init__(field=field) - self.dx = None - self.ds = None - self.n = None - self.D = None - self.S = None - self.thermal_cond = None - self.Q = None - self.T = None - self.data = [] - self.t = [] - self.show_units = False - - @property - def allowed_meshes(self): - # by default, all meshes are allowed - # override this method if that's not the case - return ["cartesian", "cylindrical", "spherical"] - - -class VolumeQuantity(DerivedQuantity): - """DerivedQuantity relative to a volume - - Args: - field (str, int): the field ("solute", 0, 1, "T", "retention") - volume (int): the volume id - - """ - - def __init__(self, field: str or int, volume: int) -> None: - super().__init__(field) - self.volume = volume - - @property - def volume(self): - return self._volume - - @volume.setter - def volume(self, value): - if not isinstance(value, int) or isinstance(value, bool): - raise TypeError("volume should be an int") - - self._volume = value - - -class SurfaceQuantity(DerivedQuantity): - """DerivedQuantity relative to a surface - - Args: - field (str, int): the field ("solute", 0, 1, "T", "retention") - surface (int): the surface id - - """ - - def __init__(self, field: str or int, surface: int) -> None: - - super().__init__(field) - self.surface = surface - - @property - def surface(self): - return self._surface - - @surface.setter - def surface(self, value): - if not isinstance(value, int) or isinstance(value, bool): - raise TypeError("surface should be an int") - self._surface = value diff --git a/exports/derived_quantities/hydrogen_flux.py b/exports/derived_quantities/hydrogen_flux.py deleted file mode 100644 index a17867580..000000000 --- a/exports/derived_quantities/hydrogen_flux.py +++ /dev/null @@ -1,26 +0,0 @@ -from festim import SurfaceFlux - - -class HydrogenFlux(SurfaceFlux): - """ - Computes the surface flux of hydrogen at a given surface - - Args: - surface (int): the surface id - - Attributes: - field (str): the hydrogen solute field - surface (int): the surface id - title (str): the title of the derived quantity - show_units (bool): show the units in the title in the derived quantities - file - function (dolfin.function.function.Function): the solution function of - the hydrogen solute field - - .. note:: - units are in H/m2/s in 1D, H/m/s in 2D and H/s in 3D domains - - """ - - def __init__(self, surface) -> None: - super().__init__(field="solute", surface=surface) diff --git a/exports/derived_quantities/maximum_surface.py b/exports/derived_quantities/maximum_surface.py deleted file mode 100644 index 4601e2faf..000000000 --- a/exports/derived_quantities/maximum_surface.py +++ /dev/null @@ -1,57 +0,0 @@ -from festim import SurfaceQuantity -import fenics as f -import numpy as np - - -class MaximumSurface(SurfaceQuantity): - """ - Computes the maximum value of a field on a given surface - - Args: - field (str, int): the field ("solute", 0, 1, "T", "retention") - surface (int): the surface id - - Attributes: - field (str, int): the field ("solute", 0, 1, "T", "retention") - surface (int): the surface id - title (str): the title of the derived quantity - show_units (bool): show the units in the title in the derived quantities - file - function (dolfin.function.function.Function): the solution function of - the field - - .. note:: - Units are in H/m3 for hydrogen concentration and K for temperature - - """ - - def __init__(self, field, surface) -> None: - super().__init__(field=field, surface=surface) - - @property - def title(self): - quantity_title = f"Maximum {self.field} surface {self.surface}" - if self.show_units: - if self.field == "T": - return quantity_title + " (K)" - else: - return quantity_title + " (H m-3)" - else: - return quantity_title - - def compute(self, surface_markers): - """Maximum of f over subdomains facets marked with self.surface""" - V = self.function.function_space() - - dm = V.dofmap() - - subd_dofs = np.unique( - np.hstack( - [ - dm.cell_dofs(c.index()) - for c in f.SubsetIterator(surface_markers, self.surface) - ] - ) - ) - - return np.max(self.function.vector().get_local()[subd_dofs]) diff --git a/exports/derived_quantities/maximum_volume.py b/exports/derived_quantities/maximum_volume.py deleted file mode 100644 index cabf82ab5..000000000 --- a/exports/derived_quantities/maximum_volume.py +++ /dev/null @@ -1,57 +0,0 @@ -from festim import VolumeQuantity -import fenics as f -import numpy as np - - -class MaximumVolume(VolumeQuantity): - """ - Computes the maximum value of a field in a given volume - - Args: - field (str, int): the field ("solute", 0, 1, "T", "retention") - volume (int): the volume id - - Attributes: - field (str, int): the field ("solute", 0, 1, "T", "retention") - volume (int): the volume id - title (str): the title of the derived quantity - show_units (bool): show the units in the title in the derived quantities - file - function (dolfin.function.function.Function): the solution function for - the field - - .. note:: - Units are in H/m3 for hydrogen concentration and K for temperature - - """ - - def __init__(self, field, volume) -> None: - super().__init__(field=field, volume=volume) - - @property - def title(self): - quantity_title = f"Maximum {self.field} volume {self.volume}" - if self.show_units: - if self.field == "T": - return quantity_title + " (K)" - else: - return quantity_title + " (H m-3)" - else: - return quantity_title - - def compute(self, volume_markers): - """Minimum of f over subdomains cells marked with self.volume""" - V = self.function.function_space() - - dm = V.dofmap() - - subd_dofs = np.unique( - np.hstack( - [ - dm.cell_dofs(c.index()) - for c in f.SubsetIterator(volume_markers, self.volume) - ] - ) - ) - - return np.max(self.function.vector().get_local()[subd_dofs]) diff --git a/exports/derived_quantities/minimum_surface.py b/exports/derived_quantities/minimum_surface.py deleted file mode 100644 index 55b576452..000000000 --- a/exports/derived_quantities/minimum_surface.py +++ /dev/null @@ -1,57 +0,0 @@ -from festim import SurfaceQuantity -import fenics as f -import numpy as np - - -class MinimumSurface(SurfaceQuantity): - """ - Computes the minimum value of a field on a given surface - - Args: - field (str, int): the field ("solute", 0, 1, "T", "retention") - surface (int): the surface id - - Attributes: - field (str, int): the field ("solute", 0, 1, "T", "retention") - surface (int): the surface id - title (str): the title of the derived quantity - show_units (bool): show the units in the title in the derived quantities - file - function (dolfin.function.function.Function): the solution function of - the field - - .. note:: - Units are in H/m3 for hydrogen concentration and K for temperature - - """ - - def __init__(self, field, surface) -> None: - super().__init__(field=field, surface=surface) - - @property - def title(self): - quantity_title = f"Minimum {self.field} surface {self.surface}" - if self.show_units: - if self.field == "T": - return quantity_title + " (K)" - else: - return quantity_title + " (H m-3)" - else: - return quantity_title - - def compute(self, surface_markers): - """Minimum of f over subdomains facets marked with self.surface""" - V = self.function.function_space() - - dm = V.dofmap() - - subd_dofs = np.unique( - np.hstack( - [ - dm.cell_dofs(c.index()) - for c in f.SubsetIterator(surface_markers, self.surface) - ] - ) - ) - - return np.min(self.function.vector().get_local()[subd_dofs]) diff --git a/exports/derived_quantities/minimum_volume.py b/exports/derived_quantities/minimum_volume.py deleted file mode 100644 index 47025118f..000000000 --- a/exports/derived_quantities/minimum_volume.py +++ /dev/null @@ -1,57 +0,0 @@ -from festim import VolumeQuantity -import fenics as f -import numpy as np - - -class MinimumVolume(VolumeQuantity): - """ - Computes the minimum value of a field in a given volume - - Args: - field (str, int): the field ("solute", 0, 1, "T", "retention") - surface (int): the surface id - - Attributes: - field (str, int): the field ("solute", 0, 1, "T", "retention") - surface (int): the surface id - title (str): the title of the derived quantity - show_units (bool): show the units in the title in the derived quantities - file - function (dolfin.function.function.Function): the solution function for - the field - - .. note:: - Units are in H/m3 for hydrogen concentration and K for temperature - - """ - - def __init__(self, field, volume) -> None: - super().__init__(field=field, volume=volume) - - @property - def title(self): - quantity_title = f"Minimum {self.field} volume {self.volume}" - if self.show_units: - if self.field == "T": - return quantity_title + " (K)" - else: - return quantity_title + " (H m-3)" - else: - return quantity_title - - def compute(self, volume_markers): - """Minimum of f over subdomains cells marked with self.volume""" - V = self.function.function_space() - - dm = V.dofmap() - - subd_dofs = np.unique( - np.hstack( - [ - dm.cell_dofs(c.index()) - for c in f.SubsetIterator(volume_markers, self.volume) - ] - ) - ) - - return np.min(self.function.vector().get_local()[subd_dofs]) diff --git a/exports/derived_quantities/point_value.py b/exports/derived_quantities/point_value.py deleted file mode 100644 index aec7fd39d..000000000 --- a/exports/derived_quantities/point_value.py +++ /dev/null @@ -1,46 +0,0 @@ -from festim import DerivedQuantity - - -class PointValue(DerivedQuantity): - """ - Computes the value of a field at a given point - - Args: - field (str, int): the field ("solute", 0, 1, "T", "retention") - x (int, float, tuple, list): the point coordinates - - Attributes: - field (str, int): the field ("solute", 0, 1, "T", "retention") - x (int, float, tuple, list): the point coordinates - title (str): the title of the derived quantity - show_units (bool): show the units in the title in the derived quantities - file - function (dolfin.function.function.Function): the solution function of - the field - - .. note:: - Units are in H/m3 for hydrogen concentration and K for temperature - - """ - - def __init__(self, field: str or int, x: int or float or tuple or list) -> None: - super().__init__(field=field) - # make sure x is an iterable - if not hasattr(x, "__iter__"): - x = [x] - self.x = x - - @property - def title(self): - quantity_title = f"{self.field} value at {self.x}" - if self.show_units: - if self.field == "T": - return quantity_title + " (K)" - else: - return quantity_title + " (H m-3)" - else: - return quantity_title - - def compute(self): - """The value at the point""" - return self.function(self.x) diff --git a/exports/derived_quantities/surface_flux.py b/exports/derived_quantities/surface_flux.py deleted file mode 100644 index 771d3a0e6..000000000 --- a/exports/derived_quantities/surface_flux.py +++ /dev/null @@ -1,292 +0,0 @@ -from festim import SurfaceQuantity, k_B -import fenics as f -import numpy as np - - -class SurfaceFlux(SurfaceQuantity): - """ - Computes the surface flux of a field at a given surface in cartesian coordinates - - Args: - field (str, int): the field ("solute", 0, 1, "T", "retention") - surface (int): the surface id - - Attributes: - field (str, int): the field ("solute", 0, 1, "T", "retention") - surface (int): the surface id - export_unit (str): the unit of the derived quantity in the export file - title (str): the title of the derived quantity - show_units (bool): show the units in the title in the derived quantities - file - function (dolfin.function.function.Function): the solution function of - the field - - .. note:: - Object to compute the flux J of a field u through a surface - J = integral(+prop * grad(u) . n ds) - where prop is the property of the field (D, thermal conductivity, etc) - u is the field - n is the normal vector of the surface - ds is the surface measure. - units are in H/m2/s in 1D, H/m/s in 2D and H/s in 3D domains for hydrogen - concentration and W/m2 in 1D, W/m in 2D and W in 3D domains for temperature - - """ - - def __init__(self, field, surface) -> None: - super().__init__(field=field, surface=surface) - self.soret = None - - @property - def allowed_meshes(self): - return ["cartesian"] - - @property - def export_unit(self): - # obtain domain dimension - dim = self.function.function_space().mesh().topology().dim() - # return unit depending on field and dimension of domain - if self.field == "T": - return f"W m{dim-3}".replace(" m0", "") - else: - return f"H m{dim-3} s-1".replace(" m0", "") - - @property - def title(self): - quantity_title = f"Flux surface {self.surface}: {self.field}" - if self.show_units: - if self.field == "T": - quantity_title = f"Heat flux surface {self.surface}" - else: - quantity_title = f"{self.field} flux surface {self.surface}" - - return quantity_title + f" ({self.export_unit})" - - else: - return quantity_title - - @property - def prop(self): - field_to_prop = { - "0": self.D, - "solute": self.D, - 0: self.D, - "T": self.thermal_cond, - } - return field_to_prop[self.field] - - def compute(self): - flux = f.assemble( - self.prop * f.dot(f.grad(self.function), self.n) * self.ds(self.surface) - ) - if self.soret and self.field in [0, "0", "solute"]: - flux += f.assemble( - self.prop - * self.function - * self.Q - / (k_B * self.T**2) - * f.dot(f.grad(self.T), self.n) - * self.ds(self.surface) - ) - return flux - - -class SurfaceFluxCylindrical(SurfaceFlux): - """ - Object to compute the flux J of a field u through a surface - J = integral(-prop * grad(u) . n ds) - where prop is the property of the field (D, thermal conductivity, etc) - u is the field - n is the normal vector of the surface - ds is the surface measure in cylindrical coordinates. - ds = r dr dtheta or ds = r dz dtheta - - .. note:: - For particle fluxes J is given in H/s, for heat fluxes J is given in W - - Args: - field (str, int): the field ("solute", 0, 1, "T", "retention") - surface (int): the surface id - azimuth_range (tuple, optional): Range of the azimuthal angle - (theta) needs to be between 0 and 2 pi. Defaults to (0, 2 * np.pi). - """ - - def __init__(self, field, surface, azimuth_range=(0, 2 * np.pi)) -> None: - super().__init__(field=field, surface=surface) - self.r = None - self.azimuth_range = azimuth_range - - @property - def export_unit(self): - # obtain domain dimension - dim = self.function.function_space().mesh().topology().dim() - # return unit depending on field and dimension of domain - if self.field == "T": - return f"W m{dim-2}".replace(" m0", "") - else: - return f"H m{dim-2} s-1".replace(" m0", "") - - @property - def allowed_meshes(self): - return ["cylindrical"] - - @property - def title(self): - if self.field == "T": - quantity_title = f"Heat flux surface {self.surface}" - else: - quantity_title = f"{self.field} flux surface {self.surface}" - - if self.show_units: - return quantity_title + f" ({self.export_unit})" - else: - return quantity_title - - @property - def azimuth_range(self): - return self._azimuth_range - - @azimuth_range.setter - def azimuth_range(self, value): - if value[0] < 0 or value[1] > 2 * np.pi: - raise ValueError("Azimuthal range must be between 0 and 2*pi") - self._azimuth_range = value - - def compute(self): - - if self.r is None: - mesh = ( - self.function.function_space().mesh() - ) # get the mesh from the function - rthetaz = f.SpatialCoordinate(mesh) # get the coordinates from the mesh - self.r = rthetaz[0] # only care about r here - - # dS_z = r dr dtheta , assuming axisymmetry dS_z = theta r dr - # dS_r = r dz dtheta , assuming axisymmetry dS_r = theta r dz - # in both cases the expression with self.ds is the same - - flux = f.assemble( - self.prop - * self.r - * f.dot(f.grad(self.function), self.n) - * self.ds(self.surface) - ) - if self.soret and self.field in [0, "0", "solute"]: - flux += f.assemble( - self.prop - * self.r - * self.function - * self.Q - / (k_B * self.T**2) - * f.dot(f.grad(self.T), self.n) - * self.ds(self.surface) - ) - flux *= self.azimuth_range[1] - self.azimuth_range[0] - return flux - - -class SurfaceFluxSpherical(SurfaceFlux): - """ - Object to compute the flux J of a field u through a surface - J = integral(-prop * grad(u) . n ds) - where prop is the property of the field (D, thermal conductivity, etc) - u is the field - n is the normal vector of the surface - ds is the surface measure in spherical coordinates. - ds = r^2 sin(theta) dtheta dphi - - .. note:: - For particle fluxes J is given in H/s, for heat fluxes J is given in W - - Args: - field (str, int): the field ("solute", 0, 1, "T", "retention") - surface (int): the surface id - azimuth_range (tuple, optional): Range of the azimuthal angle - (phi) needs to be between 0 and pi. Defaults to (0, np.pi). - polar_range (tuple, optional): Range of the polar angle - (theta) needs to be between - pi and pi. Defaults to (-np.pi, np.pi). - """ - - def __init__( - self, field, surface, azimuth_range=(0, np.pi), polar_range=(-np.pi, np.pi) - ) -> None: - super().__init__(field=field, surface=surface) - self.r = None - self.polar_range = polar_range - self.azimuth_range = azimuth_range - - @property - def export_unit(self): - if self.field == "T": - return f"W" - else: - return f"H s-1" - - @property - def allowed_meshes(self): - return ["spherical"] - - @property - def title(self): - if self.field == "T": - quantity_title = f"Heat flux surface {self.surface}" - else: - quantity_title = f"{self.field} flux surface {self.surface}" - - if self.show_units: - return quantity_title + f" ({self.export_unit})" - else: - return quantity_title - - @property - def polar_range(self): - return self._polar_range - - @polar_range.setter - def polar_range(self, value): - if value[0] < -np.pi or value[1] > np.pi: - raise ValueError("Polar range must be between - pi and pi") - self._polar_range = value - - @property - def azimuth_range(self): - return self._azimuth_range - - @azimuth_range.setter - def azimuth_range(self, value): - if value[0] < 0 or value[1] > np.pi: - raise ValueError("Azimuthal range must be between 0 and pi") - self._azimuth_range = value - - def compute(self): - - if self.r is None: - mesh = ( - self.function.function_space().mesh() - ) # get the mesh from the function - rthetaphi = f.SpatialCoordinate(mesh) # get the coordinates from the mesh - self.r = rthetaphi[0] # only care about r here - - # dS_r = r^2 sin(theta) dtheta dphi - # integral(f dS_r) = integral(f r^2 sin(theta) dtheta dphi) - # = (phi2 - phi1) * (-cos(theta2) + cos(theta1)) * f r^2 - flux = f.assemble( - self.prop - * self.r**2 - * f.dot(f.grad(self.function), self.n) - * self.ds(self.surface) - ) - if self.soret and self.field in [0, "0", "solute"]: - flux += f.assemble( - self.prop - * self.r**2 - * self.function - * self.Q - / (k_B * self.T**2) - * f.dot(f.grad(self.T), self.n) - * self.ds(self.surface) - ) - flux *= (self.polar_range[1] - self.polar_range[0]) * ( - -np.cos(self.azimuth_range[1]) + np.cos(self.azimuth_range[0]) - ) - return flux diff --git a/exports/derived_quantities/thermal_flux.py b/exports/derived_quantities/thermal_flux.py deleted file mode 100644 index 802ffd23f..000000000 --- a/exports/derived_quantities/thermal_flux.py +++ /dev/null @@ -1,26 +0,0 @@ -from festim import SurfaceFlux - - -class ThermalFlux(SurfaceFlux): - """ - Computes the surface flux of heat at a given surface - - Args: - surface (int): the surface id - - Attributes: - surface (int): the surface id - field (str): the temperature field - title (str): the title of the derived quantity - show_units (bool): show the units in the title in the derived quantities - file - function (dolfin.function.function.Function): the solution function of - the temperature field - - .. note:: - units are in W/m2 in 1D, W/m in 2D and W in 3D domains - - """ - - def __init__(self, surface) -> None: - super().__init__(field="T", surface=surface) diff --git a/exports/derived_quantities/total_surface.py b/exports/derived_quantities/total_surface.py deleted file mode 100644 index d9870ef42..000000000 --- a/exports/derived_quantities/total_surface.py +++ /dev/null @@ -1,60 +0,0 @@ -from festim import SurfaceQuantity -import fenics as f - - -class TotalSurface(SurfaceQuantity): - """ - Computes the total value of a field on a given surface - int(f ds) - - Args: - field (str, int): the field ("solute", 0, 1, "T", "retention") - surface (int): the surface id - - Attributes: - field (str, int): the field ("solute", 0, 1, "T", "retention") - surface (int): the surface id - export_unit (str): the unit of the derived quantity for exporting - title (str): the title of the derived quantity - show_units (bool): show the units in the title in the derived quantities - file - function (dolfin.function.function.Function): the solution function of - the hydrogen solute field - - .. note:: - units are in H/m2 in 1D, H/m in 2D and H in 3D domains for hydrogen - concentration and K in 1D, K m in 2D and K m2 in 3D domains for temperature - - """ - - def __init__(self, field, surface) -> None: - super().__init__(field=field, surface=surface) - - @property - def allowed_meshes(self): - return ["cartesian"] - - @property - def export_unit(self): - # obtain domain dimension - try: - dim = self.function.function_space().mesh().topology().dim() - except AttributeError: - dim = self.dx._domain._topological_dimension - # TODO we could simply do that all the time - # return unit depending on field and dimension of domain - if self.field == "T": - return f"K m{dim-1}".replace(" m0", "").replace(" m1", " m") - else: - return f"H m{dim-3}".replace(" m0", "") - - @property - def title(self): - quantity_title = f"Total {self.field} surface {self.surface}" - if self.show_units: - return quantity_title + f" ({self.export_unit})" - else: - return quantity_title - - def compute(self): - return f.assemble(self.function * self.ds(self.surface)) diff --git a/exports/derived_quantities/total_volume.py b/exports/derived_quantities/total_volume.py deleted file mode 100644 index fd770dde1..000000000 --- a/exports/derived_quantities/total_volume.py +++ /dev/null @@ -1,60 +0,0 @@ -from festim import VolumeQuantity -import fenics as f - - -class TotalVolume(VolumeQuantity): - """ - Computes the total value of a field in a given volume - int(f dx) - - Args: - field (str, int): the field ("solute", 0, 1, "T", "retention") - volume (int): the volume id - - Attributes: - field (str, int): the field ("solute", 0, 1, "T", "retention") - volume (int): the volume id - export_unit (str): the unit of the derived quantity for exporting - title (str): the title of the derived quantity - show_units (bool): show the units in the title in the derived quantities - file - function (dolfin.function.function.Function): the solution function of - the hydrogen solute field - - .. note:: - units are in H/m2 in 1D, H/m in 2D and H in 3D domains for hydrogen - concentration and K m in 1D, K m2 in 2D and K m3 in 3D domains for temperature - - """ - - def __init__(self, field, volume) -> None: - super().__init__(field=field, volume=volume) - - @property - def allowed_meshes(self): - return ["cartesian"] - - @property - def export_unit(self): - # obtain domain dimension - try: - dim = self.function.function_space().mesh().topology().dim() - except AttributeError: - dim = self.dx._domain._topological_dimension - # TODO we could simply do that all the time - # return unit depending on field and dimension of domain - if self.field == "T": - return f"K m{dim}".replace("1", "") - else: - return f"H m{dim-3}".replace(" m0", "") - - @property - def title(self): - quantity_title = f"Total {self.field} volume {self.volume}" - if self.show_units: - return quantity_title + f" ({self.export_unit})" - else: - return quantity_title - - def compute(self): - return f.assemble(self.function * self.dx(self.volume)) diff --git a/exports/exports.py b/exports/exports.py deleted file mode 100644 index 45e769a62..000000000 --- a/exports/exports.py +++ /dev/null @@ -1,154 +0,0 @@ -import festim -import fenics as f -import warnings - - -class Exports(list): - """ - A list of festim.Export objects - """ - - def __init__(self, *args): - # checks that input is list - if len(args) == 0: - super().__init__() - else: - if not isinstance(*args, list): - raise TypeError("festim.Exports must be a list") - super().__init__(self._validate_export(item) for item in args[0]) - - self.t = None - self.V_DG1 = None - self.final_time = None - self.nb_iterations = 0 - - @property - def exports(self): - warnings.warn( - "The exports attribute will be deprecated in a future release, please use festim.Exports as a list instead", - DeprecationWarning, - ) - return self - - @exports.setter - def exports(self, value): - warnings.warn( - "The exports attribute will be deprecated in a future release, please use festim.Exports as a list instead", - DeprecationWarning, - ) - if isinstance(value, list): - if not all( - ( - isinstance(t, festim.Export) - or isinstance(t, festim.DerivedQuantities) - ) - for t in value - ): - raise TypeError("exports must be a list of festim.Export") - super().__init__(value) - else: - raise TypeError("exports must be a list") - - def __setitem__(self, index, item): - super().__setitem__(index, self._validate_export(item)) - - def insert(self, index, item): - super().insert(index, self._validate_export(item)) - - def append(self, item): - super().append(self._validate_export(item)) - - def extend(self, other): - if isinstance(other, type(self)): - super().extend(other) - else: - super().extend(self._validate_export(item) for item in other) - - def _validate_export(self, value): - if isinstance(value, festim.Export) or isinstance( - value, festim.DerivedQuantities - ): - return value - raise TypeError("festim.Exports must be a list of festim.Export") - - def write(self, label_to_function, dx): - """writes to file - - Args: - label_to_function (dict): dictionary of labels mapped to solutions - dx (fenics.Measure): the measure for dx - """ - for export in self: - if isinstance(export, festim.DerivedQuantities): - # compute derived quantities - if export.is_compute(self.nb_iterations): - # check if function has to be projected - for quantity in export: - if isinstance( - quantity, (festim.MaximumVolume, festim.MinimumVolume) - ): - if not isinstance( - label_to_function[quantity.field], f.Function - ): - label_to_function[quantity.field] = f.project( - label_to_function[quantity.field], self.V_DG1 - ) - if isinstance(quantity, festim.AdsorbedHydrogen): - for surf_funcs in label_to_function[quantity.field]: - if quantity.surface in surf_funcs["surfaces"]: - ind = surf_funcs["surfaces"].index(quantity.surface) - quantity.function = surf_funcs[ - "post_processing_solutions" - ][ind] - else: - quantity.function = label_to_function[quantity.field] - export.compute(self.t) - # export derived quantities - if export.is_export(self.t, self.final_time, self.nb_iterations): - export.write() - - elif isinstance(export, festim.XDMFExport): - if export.is_export(self.t, self.final_time, self.nb_iterations): - if export.field == "retention": - # if not a Function, project it onto V_DG1 - if not isinstance(label_to_function["retention"], f.Function): - label_to_function["retention"] = f.project( - label_to_function["retention"], self.V_DG1 - ) - export.function = label_to_function[export.field] - if isinstance(export, festim.TrapDensityXDMF): - export.write(self.t, dx) - else: - export.write(self.t) - export.append = True - - elif isinstance(export, festim.TXTExport): - # if not a Function, project it onto V_DG1 - if not isinstance(label_to_function[export.field], f.Function): - label_to_function[export.field] = f.project( - 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) - self.nb_iterations += 1 - - def initialise_derived_quantities(self, dx, ds, materials): - """If derived quantities in exports, creates header and adds measures - and properties - - Args: - dx (fenics.Measure): the measure for dx - ds (fenics.Measure): the measure for ds - materials (festim.Materials): the materials - """ - for export in self: - if isinstance(export, festim.DerivedQuantities): - # reset the data of the derived quantities - export.data = [] - export.t = [] - for quantity in export: - quantity.t = [] - quantity.data = [] - export.assign_measures_to_quantities(dx, ds) - export.assign_properties_to_quantities(materials) diff --git a/exports/txt_export.py b/exports/txt_export.py deleted file mode 100644 index 3c7688a99..000000000 --- a/exports/txt_export.py +++ /dev/null @@ -1,133 +0,0 @@ -import fenics as f -import numpy as np -import festim -import warnings -import os - -warnings.simplefilter("always", DeprecationWarning) - - -class TXTExport(festim.Export): - """ - Args: - field (str): the exported field ("solute", "1", "retention", - "T"...) - filename (str): the filename (must end with .txt). - 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". - """ - - def __init__(self, field, filename, times=None, header_format=".2e") -> None: - super().__init__(field=field) - if times: - self.times = sorted(times) - else: - self.times = times - self.filename = filename - self.header_format = header_format - self._first_time = True - - @property - def filename(self): - return self._filename - - @filename.setter - def filename(self, value): - if value is not None: - if not isinstance(value, str): - raise TypeError("filename must be a string") - if not value.endswith(".txt"): - raise ValueError("filename must end with .txt") - self._filename = value - - def is_it_time_to_export(self, current_time): - if self.times is None: - return True - for time in self.times: - if np.isclose(time, current_time, atol=0): - 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 write(self, current_time, steady): - # create a DG1 functionspace - V_DG1 = f.FunctionSpace(self.function.function_space().mesh(), "DG", 1) - - solution = f.project(self.function, V_DG1) - solution_column = np.transpose(solution.vector()[:]) - if self.is_it_time_to_export(current_time): - # 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) - - # if steady or it is the first time to export - # write data - # else append new column to the existing file - if steady or self._first_time: - if steady: - header = "x,t=steady" - else: - 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._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)) diff --git a/exports/xdmf_export.py b/exports/xdmf_export.py deleted file mode 100644 index 14383fa6a..000000000 --- a/exports/xdmf_export.py +++ /dev/null @@ -1,161 +0,0 @@ -import warnings -from festim import Export -import fenics as f -import numpy as np - - -field_to_label = { - "solute": "mobile_concentration", - "T": "temperature", - "retention": "retention", - "trap": "trap_i_concentration", -} - - -class XDMFExport(Export): - """ - Args: - field (str): the exported field ("solute", "1", "retention", "T"...) - label (str, optional): label of the field in the written file. - If None, an automatic label will be given. Defaults to None. - filename (str, optional): the file path, needs to end with '.xdmf'. - If None, the label will be used. Defaults to None. - mode (int, str, optional): if "last" only the last - timestep will be exported. Otherwise the number of - iterations between each export can be provided as an integer. - Defaults to 1. - checkpoint (bool, optional): If set to True, - fenics.XDMFFile.write_checkpoint will be use, else - fenics.XDMFFile.write. Defaults to True. - folder (str, optional): path of the export folder. Defaults to None. - """ - - def __init__( - self, field, label=None, filename=None, mode=1, checkpoint=True, folder=None - ) -> None: - super().__init__(field=field) - self.label = label - self.folder = folder - self.filename = filename - - self.files = None - self.define_xdmf_file() - self.mode = mode - self.checkpoint = checkpoint - if type(self.checkpoint) != bool: - raise TypeError("checkpoint must be a bool") - - self.append = False - - @property - def label(self): - return self._label - - @label.setter - def label(self, value): - if value is None: - if self.field in field_to_label.keys(): - self._label = field_to_label[self.field] - elif self.field.isdigit(): - self._label = field_to_label["trap"].replace("i", self.field, 1) - else: - self._label = value - - @property - def mode(self): - return self._mode - - @mode.setter - def mode(self, value): - accepted_values = "accepted values for mode are int and 'last'" - if not isinstance(value, (str, int)): - raise ValueError(accepted_values) - if isinstance(value, int) and value <= 0: - raise ValueError("mode must be positive") - if isinstance(value, str) and value != "last": - raise ValueError(accepted_values) - - self._mode = value - - @property - def folder(self): - return self._folder - - @folder.setter - def folder(self, value): - if value is not None: - if not isinstance(value, str): - raise TypeError("folder must be a string") - - self._folder = value - - @property - def filename(self): - return self._filename - - @filename.setter - def filename(self, value): - if value is not None: - if not isinstance(value, str): - raise TypeError("filename must be a string") - if not value.endswith(".xdmf"): - raise ValueError("filename must end with .xdmf") - self._filename = value - else: - self._filename = "{}.xdmf".format(self.label) - - def define_xdmf_file(self): - """Creates the file""" - if self.folder is None: - filename = self.filename - else: - filename = "{}/{}".format(self.folder, self.filename) - self.file = f.XDMFFile(filename) - self.file.parameters["flush_output"] = True - self.file.parameters["rewrite_function_mesh"] = False - - def write(self, t): - """Writes to file - - Args: - t (float): current time - """ - self.function.rename(self.label, "label") - - if self.checkpoint: - # warn users if checkpoint is True and 1D - dimension = self.function.function_space().mesh().topology().dim() - if dimension == 1: - msg = "in 1D, checkpoint needs to be set to False to " - msg += "visualise the XDMF file in Paraview (see issue " - msg += "https://github.com/festim-dev/festim/issues/134)" - warnings.warn(msg) - - self.file.write_checkpoint( - self.function, - self.label, - t, - f.XDMFFile.Encoding.HDF5, - append=self.append, - ) - else: - self.file.write(self.function, t) - - def is_export(self, t, final_time, nb_iterations): - """Checks if export should be exported. - - Args: - t (float): the current time - final_time (float): the final time of the simulation - nb_iterations (int): the current number of time steps - - Returns: - bool: True if export should be exported, else False - """ - if self.mode == "last" and np.isclose(t, final_time, atol=0): - return True - elif isinstance(self.mode, int): - if nb_iterations % self.mode == 0: - return True - - return False diff --git a/generic_simulation.py b/generic_simulation.py deleted file mode 100644 index f8766d1ee..000000000 --- a/generic_simulation.py +++ /dev/null @@ -1,542 +0,0 @@ -import festim -from festim.h_transport_problem import HTransportProblem -from fenics import * -import numpy as np -import sympy as sp -import warnings - - -class Simulation: - """ - Main festim class representing a festim model - - Args: - mesh (festim.Mesh, optional): The mesh of the model. Defaults to - None. - materials (festim.Materials or list or festim.Material, optional): - The model materials. Defaults to None. - sources (list of festim.Source, optional): Volumetric sources - (particle or heat sources). Defaults to []. - boundary_conditions (list of festim.BoundaryCondition, optional): - The model's boundary conditions (temperature of H - concentration). Defaults to None. - traps (festim.Traps or list or festim.Trap, optional): The model's traps. Defaults - to None. - dt (festim.Stepsize, optional): The model's stepsize. Defaults to - None. - settings (festim.Settings, optional): The model's settings. - Defaults to None. - temperature (int, float, sympy.Expr, festim.Temperature, optional): The model's - temperature. Can be an expression or a heat transfer model. - Defaults to None. - initial_conditions (list of festim.InitialCondition, optional): - The model's initial conditions (H or T). Defaults to []. - exports (festim.Exports or list or festim.Export, optional): The model's exports - (derived quantities, XDMF exports, txt exports...). Defaults - to None. - log_level (int, optional): set what kind of FEniCS messsages are - displayed. Defaults to 40. - CRITICAL = 50, errors that may lead to data corruption - ERROR = 40, errors - WARNING = 30, warnings - INFO = 20, information of general interest - PROGRESS = 16, what's happening (broadly) - TRACE = 13, what's happening (in detail) - DBG = 10 sundry - - Attributes: - log_level (int): set what kind of FEniCS messsages are - displayed. - CRITICAL = 50, errors that may lead to data corruption - ERROR = 40, errors - WARNING = 30, warnings - INFO = 20, information of general interest - PROGRESS = 16, what's happening (broadly) - TRACE = 13, what's happening (in detail) - DBG = 10 sundry - settings (festim.Settings): The model's settings. - dt (festim.Stepsize): The model's stepsize. - traps (festim.Traps): The model's traps. - materials (festim.Materials): The model materials. - boundary_conditions (list of festim.BoundaryCondition): - The model's boundary conditions (temperature of H - concentration). - initial_conditions (list of festim.InitialCondition): - The model's initial conditions (H or T). - T (festim.Temperature): The model's temperature. - exports (festim.Exports): The model's exports - (derived quantities, XDMF exports, txt exports...). - mesh (festim.Mesh): The mesh of the model. - sources (list of festim.Source): Volumetric sources - (particle or heat sources). - mobile (festim.Mobile): the mobile concentration (c_m or theta) - t (fenics.Constant): the current time of simulation - timer (fenics.timer): the elapsed time of simulation - """ - - def __init__( - self, - mesh=None, - materials=None, - sources=[], - boundary_conditions=[], - traps=None, - dt=None, - settings=None, - temperature=None, - initial_conditions=[], - exports=None, - log_level=40, - ): - self.log_level = log_level - - self.settings = settings - self.dt = dt - - self.traps = traps - self.materials = materials - - self.boundary_conditions = boundary_conditions - self.initial_conditions = initial_conditions - self.T = temperature - self.exports = exports - self.mesh = mesh - self.sources = sources - - # internal attributes - self.h_transport_problem = None - self.t = 0 # Initialising time to 0s - self.timer = None - - @property - def traps(self): - return self._traps - - @traps.setter - def traps(self, value): - if value is None: - self._traps = festim.Traps([]) - elif isinstance(value, festim.Traps): - self._traps = value - elif isinstance(value, list): - self._traps = festim.Traps(value) - elif isinstance(value, festim.Trap): - self._traps = festim.Traps([value]) - else: - raise TypeError( - "Accepted types for traps are list, festim.Traps or festim.Trap" - ) - - @property - def materials(self): - return self._materials - - @materials.setter - def materials(self, value): - if isinstance(value, festim.Materials): - self._materials = value - elif isinstance(value, list): - self._materials = festim.Materials(value) - elif isinstance(value, festim.Material): - self._materials = festim.Materials([value]) - elif value is None: - self._materials = value - else: - raise TypeError( - "accepted types for materials are list, festim.Material or festim.Materials" - ) - - @property - def exports(self): - return self._exports - - @exports.setter - def exports(self, value): - if value is None: - self._exports = festim.Exports([]) - elif isinstance(value, (festim.Export, festim.DerivedQuantities)): - self._exports = festim.Exports([value]) - elif isinstance(value, festim.Exports): - self._exports = value - elif isinstance(value, list): - self._exports = festim.Exports(value) - else: - raise TypeError( - "accepted types for exports are list, festim.DerivedQuantities, festim.Export or festim.Exports" - ) - - @property - def T(self): - return self._T - - @T.setter - def T(self, value): - if isinstance(value, festim.Temperature): - self._T = value - elif value is None: - self._T = value - elif isinstance(value, (int, float, sp.Expr)): - self._T = festim.Temperature(value) - else: - raise TypeError( - "accepted types for T attribute are int, float, sympy.Expr or festim.Temperature" - ) - - def attribute_source_terms(self): - """Assigns the source terms (in self.sources) to the correct field - (self.mobile, self.T, or traps) - """ - # reinitialise sources for concentrations and temperature - self.mobile.sources = [] - self.T.sources = [] - for t in self.traps: - t.sources = [] - - # make field_to_object dict - field_to_object = { - "solute": self.mobile, - "0": self.mobile, - 0: self.mobile, - "mobile": self.mobile, - "T": self.T, - } - for i, trap in enumerate(self.traps, 1): - field_to_object[i] = trap - field_to_object[str(i)] = trap - - # set sources - for source in self.sources: - if source.field == "T" and not isinstance( - self.T, festim.HeatTransferProblem - ): # check that there is not a source defined in T as the same time as a festim.Temperature - raise TypeError( - "Heat transfer sources can only be used with HeatTransferProblem" - ) - if isinstance(source, festim.RadioactiveDecay) and source.field == "all": - # assign source to each of the unique festim.Concentration - # objects in field_to_object - for obj in set(field_to_object.values()): - if isinstance(obj, festim.Concentration): - obj.sources.append(source) - else: - field_to_object[source.field].sources.append(source) - - def check_boundary_conditions(self): - """Runs a series of checks on the BCs and raise errors accordingly""" - - valid_fields = ( - ["T", 0, "0"] # temperature and mobile concentration - + [str(i + 1) for i, _ in enumerate(self.traps)] - + [i + 1 for i, _ in enumerate(self.traps)] - ) - - # collect all DirichletBCs and SurfaceKinetics objects - dc_sk_bcs = [ - bc - for bc in self.boundary_conditions - if isinstance(bc, (festim.DirichletBC, festim.SurfaceKinetics)) - ] - - for bc in self.boundary_conditions: - if bc.field not in valid_fields: - raise ValueError(f"{bc.field} is not a valid field for BC") - - # check SurfaceKinetics in 1D simulations - if ( - isinstance(bc, festim.SurfaceKinetics) - and self.mesh.mesh.topology().dim() != 1 - ): - raise ValueError("SurfaceKinetics can only be used in 1D simulations") - - # check that there is not a Temperature defined at the same time as a boundary condition in T - if bc.field == "T" and not isinstance(self.T, festim.HeatTransferProblem): - raise TypeError( - "Heat transfer boundary conditions can only be used with HeatTransferProblem" - ) - - # checks that DirichletBC or SurfaceKinetics is not set with another bc on the same surface - # iterate through all BCs - for dc_sk_bc in dc_sk_bcs: - if ( - bc == dc_sk_bc or bc.field != dc_sk_bc.field - ): # skip if the same BC or different fields - continue - - # check if BCs share the same surfaces using the set().isdisjoint() method - # that returns True if the first set has no elements in common with other containers - if not set(bc.surfaces).isdisjoint(dc_sk_bc.surfaces): - # convert lists of surfaces to sets and obtain their intersection - intersection = set(bc.surfaces) & set(dc_sk_bc.surfaces) - - # check the bc type for the export message - bc_type = ( - "DirichletBC" - if isinstance(dc_sk_bc, festim.DirichletBC) - else "SurfaceKinetics" - ) - msg = f"{bc_type} is simultaneously set with another boundary condition " - msg += f"on surfaces {intersection} for field {dc_sk_bc.field}" - raise ValueError(msg) - - def attribute_boundary_conditions(self): - """Assigns boundary_conditions to mobile and T""" - self.T.boundary_conditions = [] - self.h_transport_problem.boundary_conditions = [] - self.check_boundary_conditions() - - for bc in self.boundary_conditions: - if bc.field == "T": - self.T.boundary_conditions.append(bc) - else: - self.h_transport_problem.boundary_conditions.append(bc) - - def initialise(self): - """Initialise the model. Defines markers, create the suitable function - spaces, the functions, the variational forms... - """ - set_log_level(self.log_level) - - self.t = 0 # reinitialise t to zero - - if self.settings.chemical_pot: - self.mobile = festim.Theta() - else: - self.mobile = festim.Mobile() - # check that dt attribute is None if the sim is steady state - if not self.settings.transient and self.dt is not None: - raise AttributeError("dt must be None in steady state simulations") - if self.settings.transient and self.settings.final_time is None: - raise AttributeError( - "final_time argument must be provided to settings in transient simulations" - ) - if self.settings.transient and self.dt is None: - raise AttributeError("dt must be provided in transient simulations") - if not self.T: - raise AttributeError("Temperature is not defined") - - # initialise dt - if self.settings.transient: - self.dt.initialise_value() - - self.h_transport_problem = HTransportProblem( - self.mobile, self.traps, self.T, self.settings, self.initial_conditions - ) - self.attribute_source_terms() - self.attribute_boundary_conditions() - - if isinstance(self.mesh, festim.Mesh1D): - self.mesh.define_measures(self.materials) - else: - self.mesh.define_measures() - - # needed to avoid hanging behaviour in parrallel see #498 - self.mesh.mesh.bounding_box_tree() - - self.V_DG1 = FunctionSpace(self.mesh.mesh, "DG", 1) - self.exports.V_DG1 = self.V_DG1 - - # Define temperature - if isinstance(self.T, festim.HeatTransferProblem): - self.T.create_functions(self.materials, self.mesh, self.dt) - elif isinstance(self.T, festim.Temperature): - self.T.create_functions(self.mesh) - - # Create functions for properties - self.materials.check_materials( - self.T, derived_quantities=[] - ) # FIXME derived quantities shouldn't be [] - self.materials.create_properties(self.mesh.volume_markers, self.T.T) - self.materials.create_solubility_law_markers(self.mesh) - - # if the temperature is not time-dependent, solubility can be projected - if self.settings.chemical_pot: - # TODO this could be moved to Materials.create_properties() - if self.T.is_steady_state(): - # self.materials.S = project(self.materials.S, self.V_DG1) - self.materials.solubility_as_function(self.mesh, self.T.T) - - 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: - if self.mesh.type not in q.allowed_meshes: - warnings.warn( - f"{type(q)} may not work as intended for {self.mesh.type} meshes" - ) - - if isinstance(q, festim.AdsorbedHydrogen): - # check that festim.AdsorbedHydrogen is defined together with - # festim.SurfaceKinetics on the same surface - surf_kin_present = any( - q.surface in bc.surfaces - for bc in self.boundary_conditions - if isinstance(bc, festim.SurfaceKinetics) - ) - - if not surf_kin_present: - raise AttributeError( - 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: - if isinstance(q, festim.SurfaceFlux): - q.soret = self.settings.soret - q.T = self.T.T - - def run(self, completion_tone=False): - """Runs the model. - - Args: - completion_tone (bool, optional): If True, a native os alert - tone will alert user upon completion of current run. Defaults - to False. - Returns: - dict: output containing solutions, mesh, derived quantities - """ - self.timer = Timer() # start timer - - if self.settings.transient: - self.run_transient() - else: - self.run_steady() - - self.timer.stop() - - # End - if completion_tone: - print("\007") - - def run_transient(self): - # add final_time to Exports - self.exports.final_time = self.settings.final_time - - # compute Jacobian before iterating if required - if not self.settings.update_jacobian: - self.h_transport_problem.compute_jacobian() - - # Time-stepping - print("Time stepping...") - while self.t < self.settings.final_time and not np.isclose( - self.t, self.settings.final_time, atol=0 - ): - self.iterate() - - def run_steady(self): - # Solve steady state - print("Solving steady state problem...") - - nb_iterations, converged = self.h_transport_problem.solve_once() - - # Post processing - self.run_post_processing() - elapsed_time = round(self.timer.elapsed()[0], 1) - - # print final message - if converged: - msg = "Solved problem in {:.2f} s".format(elapsed_time) - print(msg) - else: - msg = "The solver diverged in " - msg += "{:.0f} iteration(s) ({:.2f} s)".format(nb_iterations, elapsed_time) - raise ValueError(msg) - - def iterate(self): - """Advance the model by one iteration""" - # Update current time - self.t += float(self.dt.value) - # update temperature - self.T.update(self.t) - # update H problem - self.h_transport_problem.update(self.t, self.dt) - - # Display time - self.display_time() - - # Post processing - self.run_post_processing() - - # avoid t > final_time - next_time = self.t + float(self.dt.value) - if next_time > self.settings.final_time: - self.dt.value.assign(self.settings.final_time - self.t) - - def display_time(self): - """Displays the current time""" - simulation_percentage = round(self.t / self.settings.final_time * 100, 2) - elapsed_time = round(self.timer.elapsed()[0], 1) - msg = "{:.1f} % ".format(simulation_percentage) - msg += "{:.1e} s".format(self.t) - msg += " Elapsed time so far: {:.1f} s".format(elapsed_time) - if ( - not np.isclose(self.t, self.settings.final_time, atol=0) - and self.log_level == 40 - ): - print(msg, end="\r") - else: - print(msg) - - def run_post_processing(self): - """Create post processing functions and compute/write the exports""" - self.update_post_processing_solutions() - - self.exports.t = self.t - self.exports.write(self.label_to_function, self.mesh.dx) - - def update_post_processing_solutions(self): - """Creates the post-processing functions by splitting self.u. Projects - the function on a suitable functionspace if needed. - - Returns: - dict: a mapping of the field ("solute", "T", "retention") to its - post_processsing_solution - """ - self.h_transport_problem.update_post_processing_solutions(self.exports) - - label_to_function = { - "solute": self.mobile.post_processing_solution, - "0": self.mobile.post_processing_solution, - 0: self.mobile.post_processing_solution, - "T": self.T.T, - "retention": sum( - [self.mobile.post_processing_solution] - + [trap.post_processing_solution for trap in self.traps] - ), - # dictionary {"post_processing_solutions": bc.post_processing_solutions, "surfaces": bc.surfaces} - # for each SurfaceKinetics boundary condition - "adsorbed": [ - { - "post_processing_solutions": bc.post_processing_solutions, - "surfaces": bc.surfaces, - } - for bc in self.boundary_conditions - if isinstance(bc, festim.SurfaceKinetics) - ], - } - for trap in self.traps: - label_to_function[trap.id] = trap.post_processing_solution - label_to_function[str(trap.id)] = trap.post_processing_solution - - self.label_to_function = label_to_function diff --git a/h_transport_problem.py b/h_transport_problem.py deleted file mode 100644 index 7d4c4859f..000000000 --- a/h_transport_problem.py +++ /dev/null @@ -1,388 +0,0 @@ -from fenics import * -import festim - - -class HTransportProblem: - """Hydrogen Transport Problem. - Used internally in festim.Simulation - - Args: - mobile (festim.Mobile): the mobile concentration - traps (festim.Traps): the traps - T (festim.Temperature): the temperature - settings (festim.Settings): the problem settings - initial_conditions (list of festim.initial_conditions): the - initial conditions of the h transport problem - - Attributes: - expressions (list): contains time-dependent fenics.Expressions - J (ufl.Form): the jacobian of the variational problem - V (fenics.FunctionSpace): the vector-function space for concentrations - u (fenics.Function): the vector holding the concentrations (c_m, ct1, - ct2, ...) - v (fenics.TestFunction): the test function - u_n (fenics.Function): the "previous" function - newton_solver (fenics.NewtonSolver): Newton solver for solving the nonlinear problem - bcs (list): list of fenics.DirichletBC for H transport - """ - - def __init__(self, mobile, traps, T, settings, initial_conditions) -> None: - self.mobile = mobile - self.traps = traps - self.T = T - self.settings = settings - self.initial_conditions = initial_conditions - - self.J = None - self.u = None - self.v = None - self.u_n = None - self.newton_solver = None - - self.boundary_conditions = [] - self.bcs = None - self.V = None - self.V_CG1 = None - self.expressions = [] - - @property - def newton_solver(self): - return self._newton_solver - - @newton_solver.setter - def newton_solver(self, value): - if value is None: - self._newton_solver = value - elif isinstance(value, NewtonSolver): - if self._newton_solver: - print("Settings for the Newton solver will be overwritten") - self._newton_solver = value - else: - raise TypeError("accepted type for newton_solver is fenics.NewtonSolver") - - @property - def _all_surf_kinetics(self): - return [ - bc - for bc in self.boundary_conditions - if isinstance(bc, festim.SurfaceKinetics) - ] - - def initialise(self, mesh, materials, dt=None): - """Assigns BCs, create suitable function space, initialise - concentration fields, define variational problem - - Args: - mesh (festim.Mesh): the mesh - materials (festim.Materials): the materials - dt (festim.Stepsize, optional): the stepsize, only needed if - self.settings.transient is True. Defaults to None. - """ - if self.settings.chemical_pot: - self.mobile.S = materials.S - self.mobile.materials = materials - self.mobile.volume_markers = mesh.volume_markers - self.mobile.T = self.T - self.attribute_flux_boundary_conditions() - - self.traps.assign_traps_ids() - - # Define functions - self.define_function_space(mesh) - self.initialise_concentrations() - self.traps.make_traps_materials(materials) - self.traps.initialise_extrinsic_traps(self.V_CG1) - - # Define variational problem H transport - # if chemical pot create form to convert theta to concentration - if self.settings.chemical_pot: - self.mobile.create_form_post_processing(self.V_DG1, materials, mesh.dx) - - self.define_variational_problem(materials, mesh, dt) - self.define_newton_solver() - - # Boundary conditions - print("Defining boundary conditions") - self.create_dirichlet_bcs(materials, mesh) - if self.settings.transient: - self.traps.define_variational_problem_extrinsic_traps(mesh.dx, dt, self.T) - self.traps.define_newton_solver_extrinsic_traps() - - def define_function_space(self, mesh): - """Creates a suitable function space for H transport problem - - Args: - mesh (festim.Mesh): the mesh - """ - order_trap = 1 - element_solute, order_solute = "CG", 1 - - # function space for H concentrations - nb_traps = len(self.traps) - - # the number of surfaces where SurfaceKinetics is used - nb_adsorbed = sum([len(bc.surfaces) for bc in self._all_surf_kinetics]) - - if nb_traps == 0 and nb_adsorbed == 0: - V = FunctionSpace(mesh.mesh, element_solute, order_solute) - else: - solute = FiniteElement(element_solute, mesh.mesh.ufl_cell(), order_solute) - traps = FiniteElement( - self.settings.traps_element_type, mesh.mesh.ufl_cell(), order_trap - ) - adsorbed = FiniteElement("R", mesh.mesh.ufl_cell(), 0) - element = [solute] + [traps] * nb_traps + [adsorbed] * nb_adsorbed - V = FunctionSpace(mesh.mesh, MixedElement(element)) - - self.V = V - self.V_CG1 = FunctionSpace(mesh.mesh, "CG", 1) - self.V_DG1 = FunctionSpace(mesh.mesh, "DG", 1) - - def initialise_concentrations(self): - """Creates the main fenics.Function (holding all the concentrations), - eventually split it and assign it to Trap and Mobile. - Then initialise self.u_n based on self.initial_conditions - - Args: - materials (festim.Materials): the materials - """ - # TODO rename u and u_n to c and c_n - self.u = Function(self.V, name="c") # Function for concentrations - self.v = TestFunction(self.V) # TestFunction for concentrations - self.u_n = Function(self.V, name="c_n") - - if self.V.num_sub_spaces() == 0: - self.mobile.solution = self.u - self.mobile.previous_solution = self.u_n - self.mobile.test_function = self.v - else: - conc_list = [self.mobile] - if self.traps: - conc_list += [*self.traps] - if len(self._all_surf_kinetics) > 0: - conc_list += self._all_surf_kinetics - - index = 0 - for concentration in conc_list: - if isinstance(concentration, festim.SurfaceKinetics): - # iterate through each surface of each SurfaceKinetics - for i in range(len(concentration.surfaces)): - concentration.solutions[i] = self.u.sub(index) - concentration.previous_solutions[i] = self.u_n.sub(index) - concentration.test_functions[i] = list(split(self.v))[index] - index += 1 - else: - concentration.solution = self.u.sub(index) - concentration.previous_solution = self.u_n.sub(index) - concentration.test_function = list(split(self.v))[index] - index += 1 - - print("Defining initial values") - field_to_component = { - "solute": 0, - "0": 0, - 0: 0, - } - for i, trap in enumerate(self.traps, 1): - field_to_component[trap.id] = i - field_to_component[str(trap.id)] = i - # TODO refactore this, attach the initial conditions to the objects directly - for ini in self.initial_conditions: - value = ini.value - component = field_to_component[ini.field] - - if self.V.num_sub_spaces() == 0: - functionspace = self.V - else: - functionspace = self.V.sub(component).collapse() - - if component == 0: - self.mobile.initialise( - functionspace, value, label=ini.label, time_step=ini.time_step - ) - else: - trap = self.traps.get_trap(component) - trap.initialise( - functionspace, value, label=ini.label, time_step=ini.time_step - ) - - # assign initial condition for SurfaceKinetics BC - # iterate through each surface of each SurfaceKinetics - index = len(self.traps) + 1 - for bc in self._all_surf_kinetics: - for i in range(len(bc.previous_solutions)): - functionspace = self.V.sub(index).collapse() - comp = interpolate(Constant(bc.initial_condition), functionspace) - assign(bc.previous_solutions[i], comp) - index += 1 - - # initial guess needs to be non zero if chemical pot - if self.settings.chemical_pot: - if self.V.num_sub_spaces() == 0: - functionspace = self.V - else: - functionspace = self.V.sub(0).collapse() - initial_guess = project( - self.mobile.previous_solution + Constant(DOLFIN_EPS), functionspace - ) - self.mobile.solution.assign(initial_guess) - # this is needed to correctly create the formulation - # TODO: write a test for this? - if self.V.num_sub_spaces() != 0: - index = 0 - for concentration in conc_list: - if isinstance(concentration, festim.SurfaceKinetics): - for i in range(len(concentration.surfaces)): - concentration.solutions[i] = list(split(self.u))[index] - concentration.previous_solutions[i] = list(split(self.u_n))[ - index - ] - index += 1 - else: - concentration.solution = list(split(self.u))[index] - concentration.previous_solution = list(split(self.u_n))[index] - index += 1 - - def define_variational_problem(self, materials, mesh, dt=None): - """Creates the variational problem for hydrogen transport (form, - Dirichlet boundary conditions) - - Args: - materials (festim.Materials): the materials - mesh (festim.Mesh): the mesh - dt (festim.Stepsize, optional): the stepsize, only needed if - self.settings.transient is True. Defaults to None. - """ - print("Defining variational problem") - expressions = [] - F = 0 - - # diffusion + transient terms - - self.mobile.create_form( - materials, mesh, self.T, dt, traps=self.traps, soret=self.settings.soret - ) - F += self.mobile.F - expressions += self.mobile.sub_expressions - - # Add traps - self.traps.create_forms(self.mobile, materials, self.T, mesh.dx, dt) - F += self.traps.F - expressions += self.traps.sub_expressions - self.F = F - self.expressions = expressions - - def define_newton_solver(self): - """Creates the Newton solver and sets its parameters""" - self.newton_solver = NewtonSolver(MPI.comm_world) - self.newton_solver.parameters["error_on_nonconvergence"] = False - self.newton_solver.parameters["absolute_tolerance"] = ( - self.settings.absolute_tolerance - ) - self.newton_solver.parameters["relative_tolerance"] = ( - self.settings.relative_tolerance - ) - self.newton_solver.parameters["maximum_iterations"] = ( - self.settings.maximum_iterations - ) - self.newton_solver.parameters["linear_solver"] = self.settings.linear_solver - self.newton_solver.parameters["preconditioner"] = self.settings.preconditioner - - def attribute_flux_boundary_conditions(self): - """Iterates through self.boundary_conditions, checks if it's a FluxBC - and its field is 0, and assign fluxes to self.mobile - """ - for bc in self.boundary_conditions: - if isinstance(bc, festim.FluxBC) and bc.field == 0: - self.mobile.boundary_conditions.append(bc) - - def create_dirichlet_bcs(self, materials, mesh): - """Creates fenics.DirichletBC objects for the hydrogen transport - problem and add them to self.bcs - """ - self.bcs = [] - for bc in self.boundary_conditions: - if bc.field != "T" and isinstance(bc, festim.DirichletBC): - bc.create_dirichletbc( - self.V, - self.T.T, - mesh.surface_markers, - chemical_pot=self.settings.chemical_pot, - materials=materials, - volume_markers=mesh.volume_markers, - ) - self.bcs += bc.dirichlet_bc - self.expressions += bc.sub_expressions - self.expressions.append(bc.expression) - - def compute_jacobian(self): - du = TrialFunction(self.u.function_space()) - self.J = derivative(self.F, self.u, du) - - def update(self, t, dt): - """Updates the H transport problem. - - Args: - t (float): the current time (s) - dt (festim.Stepsize): the stepsize - """ - festim.update_expressions(self.expressions, t) - - converged = False - u_ = Function(self.u.function_space()) - u_.assign(self.u) - while converged is False: - self.u.assign(u_) - nb_it, converged = self.solve_once() - if dt.adaptive_stepsize is not None or dt.milestones is not None: - dt.adapt(t, nb_it, converged) - - # Update previous solutions - self.update_previous_solutions() - - # Solve extrinsic traps formulation - self.traps.solve_extrinsic_traps() - - def solve_once(self): - """Solves non linear problem - - Returns: - int, bool: number of iterations for reaching convergence, True if - converged else False - """ - if self.J is None: # Define the Jacobian - du = TrialFunction(self.u.function_space()) - J = derivative(self.F, self.u, du) - else: - J = self.J - problem = festim.Problem(J, self.F, self.bcs) - - begin("Solving nonlinear variational problem.") # Add message to fenics logs - nb_it, converged = self.newton_solver.solve(problem, self.u.vector()) - end() - - return nb_it, converged - - def update_previous_solutions(self): - self.u_n.assign(self.u) - self.traps.update_extrinsic_traps_density() - - def update_post_processing_solutions(self, exports): - if self.u.function_space().num_sub_spaces() == 0: - res = [self.u] - else: - res = list(self.u.split()) - - for i, trap in enumerate(self.traps, 1): - trap.post_processing_solution = res[i] - - index = len(self.traps) + 1 - for bc in self._all_surf_kinetics: - for i in range(len(bc.post_processing_solutions)): - bc.post_processing_solutions[i] = res[index] - index += 1 - - if self.settings.chemical_pot: - self.mobile.post_processing_solution_to_concentration() - else: - self.mobile.post_processing_solution = res[0] diff --git a/initial_condition.py b/initial_condition.py deleted file mode 100644 index 80065f4c2..000000000 --- a/initial_condition.py +++ /dev/null @@ -1,27 +0,0 @@ -class InitialCondition: - """ - Args: - field (int, str, optional): the field - ("0", "solute", "T", "1",...). Defaults to 0. - value (float, str, optional): the value of the initial condition. - Defaults to 0. - label (str, optional): label in the XDMF file. Defaults to None. - time_step ([type], optional): [description]. Defaults to None. - - Raises: - ValueError: if XDMF and label is None - ValueError: if XDMF and time_step is None - """ - - def __init__(self, field=0, value=0.0, label=None, time_step=None) -> None: - # TODO make an inherited class InitialConditionXDMF - self.field = field - self.value = value - self.label = label - self.time_step = time_step - if type(self.value) == str: - if self.value.endswith(".xdmf"): - if self.label is None: - raise ValueError("label cannot be None") - if self.time_step is None: - raise ValueError("time_step cannot be None") diff --git a/materials/material.py b/materials/material.py deleted file mode 100644 index 5e1ee359b..000000000 --- a/materials/material.py +++ /dev/null @@ -1,83 +0,0 @@ -class Material: - """ - Args: - id (int, list): the id of the material. If a list is provided, the - properties will be applied to all the subdomains with the - corresponding ids. - D_0 (float): diffusion coefficient pre-exponential factor (m2/s) - E_D (float): diffusion coefficient activation energy (eV) - S_0 (float, optional): Solubility pre-exponential factor - (H/m3/Pa0.5). Defaults to None. - E_S (float, optional): Solubility activation energy (eV). - Defaults to None. - thermal_cond (float or callable, optional): thermal conductivity - (W/m/K). Can be a function of T. Defaults to None. - heat_capacity (float or callable, optional): heat capacity - (J/K/kg). Can be a function of T. Defaults to None. - rho (float or callable, optional): volumetric density (kg/m3). Can - be a function of T. Defaults to None. - borders (list, optional): The borders of the 1D subdomain. - Only needed in 1D with several materials. Defaults to None. - Q (float or callable, optional): heat of transport (eV). Can - be a function of T. Defaults to None. - solubility_law (str, optional): the material's solubility law. - Can be "henry" or "sievert". Defaults to "sievert". - name (str, optional): name of the material. Defaults to None. - - Example:: - - my_mat = Material( - id=1, - D_0=2e-7, - E_d=0.2, - thermal_cond=lambda T: 3 * T + 2, - heat_capacity=lambda T: 4 * T + 8, - rho=lambda T: 7 * T + 5, - Q=lambda T: -0.5 * T**2, - ) - """ - - def __init__( - self, - id, - D_0, - E_D, - S_0=None, - E_S=None, - thermal_cond=None, - heat_capacity=None, - rho=None, - borders=None, - Q=None, - solubility_law="sievert", - name=None, - ) -> None: - self.id = id - self.name = name - self.D_0 = D_0 - self.E_D = E_D - self.S_0 = S_0 - self.E_S = E_S - self.thermal_cond = thermal_cond - self.heat_capacity = heat_capacity - self.rho = rho - self.borders = borders - self.Q = Q - if solubility_law not in ["henry", "sievert"]: - raise ValueError( - "Acceptable values for solubility_law are 'henry' and 'sievert'" - ) - self.solubility_law = solubility_law - self.check_properties() - - def check_properties(self): - """Checks that if S_0 is None E_S is not None and reverse. - - Raises: - ValueError: [description] - ValueError: [description] - """ - if self.S_0 is None and self.E_S is not None: - raise ValueError("S_0 cannot be None") - if self.E_S is None and self.S_0 is not None: - raise ValueError("E_S cannot be None") diff --git a/materials/materials.py b/materials/materials.py deleted file mode 100644 index 642f6dde4..000000000 --- a/materials/materials.py +++ /dev/null @@ -1,410 +0,0 @@ -from operator import itemgetter -import warnings -import numpy as np -from festim import k_B, Material, HeatTransferProblem -import festim -import fenics as f -from typing import Union -import warnings - - -class Materials(list): - """ - A list of festim.Material objects - """ - - def __init__(self, *args): - # checks that input is list - if len(args) == 0: - super().__init__() - else: - if not isinstance(*args, list): - raise TypeError("festim.Materials must be a list") - super().__init__(self._validate_material(item) for item in args[0]) - - self.D = None - self.S = None - self.thermal_cond = None - self.heat_capacity = None - self.density = None - self.Q = None - - @property - def materials(self): - warnings.warn( - "The materials attribute will be deprecated in a future release, please use festim.Materials as a list instead", - DeprecationWarning, - ) - return self - - @materials.setter - def materials(self, value): - warnings.warn( - "The materials attribute will be deprecated in a future release, please use festim.Materials as a list instead", - DeprecationWarning, - ) - if isinstance(value, list): - if not all(isinstance(t, festim.Material) for t in value): - raise TypeError("materials must be a list of festim.Material") - super().__init__(value) - else: - raise TypeError("materials must be a list") - - def __setitem__(self, index, item): - super().__setitem__(index, self._validate_material(item)) - - def insert(self, index, item): - super().insert(index, self._validate_material(item)) - - def append(self, item): - super().append(self._validate_material(item)) - - def extend(self, other): - if isinstance(other, type(self)): - super().extend(other) - else: - super().extend(self._validate_material(item) for item in other) - - def _validate_material(self, value): - if isinstance(value, festim.Material): - return value - raise TypeError("festim.Materials must be a list of festim.Material") - - def check_borders(self, size): - """Checks that the borders of the materials match - - Args: - size (float): size of the 1D domain - - Raises: - ValueError: if the borders don't begin at zero - ValueError: if borders don't match - ValueError: if borders don't end at size - - Returns: - bool -- True if everything's alright - """ - all_borders = [] - for m in self: - if isinstance(m.borders[0], list): - for border in m.borders: - all_borders.append(border) - else: - all_borders.append(m.borders) - all_borders = sorted(all_borders, key=itemgetter(0)) - if all_borders[0][0] != 0: - raise ValueError("Borders don't begin at zero") - for i in range(0, len(all_borders) - 1): - if all_borders[i][1] != all_borders[i + 1][0]: - raise ValueError("Borders don't match to each other") - if all_borders[len(all_borders) - 1][1] != size: - raise ValueError("Borders don't match with size") - return True - - def check_materials(self, T: festim.Temperature, derived_quantities: list = []): - """Checks the materials keys - - Args: - T (festim.Temperature): the temperature - derived_quantities (list): list of festim.DerivedQuantity - objects the derived quantities. Defaults to []. - """ - - if len(self) > 0: # TODO: get rid of this... - self.check_consistency() - - self.check_for_unused_properties(T, derived_quantities) - - self.check_unique_ids() - - self.check_missing_properties(T, derived_quantities) - - def check_unique_ids(self): - # check that ids are different - mat_ids = [] - for mat in self: - if type(mat.id) is list: - mat_ids += mat.id - else: - mat_ids.append(mat.id) - - if len(mat_ids) != len(np.unique(mat_ids)): - raise ValueError("Some materials have the same id") - - def check_for_unused_properties( - self, T: festim.Temperature, derived_quantities: list - ): - """Warns users if properties will be ignored - - Args: - T (festim.Temperature): the temperature - derived_quantities (list): list of festim.DerivedQuantity - objects - """ - # TODO add a check for ignored solubility when chemical_pot is False - # warn about unused keys - transient_properties = ["rho", "heat_capacity"] - if not isinstance(T, HeatTransferProblem): - for mat in self: - for key in transient_properties: - if getattr(mat, key) is not None: - warnings.warn(key + " key will be ignored", UserWarning) - - for mat in self: - if getattr(mat, "thermal_cond") is not None: - warn = True - if isinstance(T, HeatTransferProblem): - warn = False - else: - surface_fluxes = list( - quant - for quant in derived_quantities - if isinstance(quant, festim.SurfaceFlux) - ) - - for surface_flux in surface_fluxes: - if surface_flux.field == "T": - warn = False - if warn: - warnings.warn("thermal_cond key will be ignored", UserWarning) - - def check_consistency(self): - """Checks that materials have the same attributes""" - # check the materials keys match - attributes = { - "S_0": [], - "E_S": [], - "thermal_cond": [], - "heat_capacity": [], - "rho": [], - "borders": [], - "Q": [], - } - - for attr, value in attributes.items(): - for mat in self: - value.append(getattr(mat, attr)) - if value.count(None) not in [0, len(self)]: - raise ValueError("{} is not defined for all materials".format(attr)) - - def check_missing_properties(self, T: festim.Temperature, derived_quantities: list): - """Checks if the materials miss some properties - - Args: - T (festim.Temperature): the temperature - derived_quantities (list): list of festim.DerivedQuantity objects - - Raises: - ValueError: if thermal_cond, heat_capacity or rho is None when needed - """ - if isinstance(T, HeatTransferProblem): - if self[0].thermal_cond is None: - raise ValueError("Missing thermal_cond in materials") - if T.transient: - if self[0].heat_capacity is None: - raise ValueError("Missing heat_capacity in materials") - if self[0].rho is None: - raise ValueError("Missing rho in materials") - # TODO: add check for thermal cond for thermal flux computation - - def find_material_from_id(self, mat_id): - """Returns the material from a given id - - Args: - mat_id (int): id of the wanted material - - Raises: - ValueError: if the id isn't found - - Returns: - festim.Material: the material that has the id mat_id - """ - for material in self: - mat_ids = material.id - if type(mat_ids) is not list: - mat_ids = [mat_ids] - if mat_id in mat_ids: - return material - raise ValueError("Couldn't find ID " + str(mat_id) + " in materials list") - - def find_material_from_name(self, name): - """Returns the material with the correct name - - Args: - name (str): the name of the material - - Raises: - ValueError: when no match was found - - Returns: - festim.Material: the material object - """ - for material in self: - if material.name == name: - return material - - msg = "No material with name {} was found".format(name) - raise ValueError(msg) - - def find_material(self, mat: Union[int, str, Material]): - """Returns the correct festim.Material object based on either an id, - a name - - Args: - mat (Union[int, str, Material]): the material tag - - Returns: - festim.Material: the matching material - """ - if isinstance(mat, int): - return self.find_material_from_id(mat) - elif isinstance(mat, str): - return self.find_material_from_name(mat) - elif isinstance(mat, Material): - return mat - - def find_subdomain_from_x_coordinate(self, x): - """Finds the correct subdomain at a given x coordinate - - Args: - x (float): the x coordinate - - Returns: - int: the corresponding subdomain id - """ - for material in self: - # if no borders are provided, assume only one subdomain - if material.borders is None: - return material.id - # else find the correct material - else: - if isinstance(material.borders[0], list) and len(material.borders) > 1: - list_of_borders = material.borders - else: - list_of_borders = [material.borders] - if isinstance(material.id, list): - subdomains = material.id - else: - subdomains = [material.id for _ in range(len(list_of_borders))] - - for borders, subdomain in zip(list_of_borders, subdomains): - if borders[0] <= x <= borders[1]: - return subdomain - # if no subdomain was found, return 0 - return 0 - - def create_properties(self, vm, T): - """Creates the properties fields needed for post processing - - Arguments: - vm {fenics.MeshFunction()} -- volume markers - T {fenics.Function()} -- temperature - """ - self.D = ArheniusCoeff(self, vm, T, "D_0", "E_D", degree=2) - # all materials have the same properties so only checking the first is enough - if self[0].S_0 is not None: - self.S = ArheniusCoeff(self, vm, T, "S_0", "E_S", degree=2) - if self[0].thermal_cond is not None: - self.thermal_cond = ThermalProp(self, vm, T, "thermal_cond", degree=2) - self.heat_capacity = ThermalProp(self, vm, T, "heat_capacity", degree=2) - self.density = ThermalProp(self, vm, T, "rho", degree=2) - if self[0].Q is not None: - self.Q = ThermalProp(self, vm, T, "Q", degree=2) - - def solubility_as_function(self, mesh, T): - """ - Makes solubility as a fenics.Function and stores it in S attribute - """ - V = f.FunctionSpace(mesh.mesh, "DG", 1) - S = f.Function(V, name="S") - vS = f.TestFunction(V) - dx = mesh.dx - F = 0 - for mat in self: - F += -S * vS * dx(mat.id) - F += mat.S_0 * f.exp(-mat.E_S / k_B / T) * vS * dx(mat.id) - f.solve(F == 0, S, bcs=[]) - - self.S = S - - def create_solubility_law_markers(self, mesh: festim.Mesh): - """Creates the attributes henry_marker and sievert_marker - These fenics.Function are equal to one or zero depending - on the material solubility_law - - Args: - mesh (festim.Mesh): the mesh - """ - V = f.FunctionSpace(mesh.mesh, "DG", 0) - henry = f.Function(V) - sievert = f.Function(V) - - test_function_henry = f.TestFunction(V) - test_function_sievert = f.TestFunction(V) - - # initialise formulations - F_henry = -henry * test_function_henry * mesh.dx - F_sievert = -sievert * test_function_sievert * mesh.dx - - # build the formulation depending on the - for mat in self: - # make sure mat_ids is a list - mat_ids = mat.id - if not isinstance(mat.id, list): - mat_ids = [mat.id] - - for mat_id in mat_ids: # iterate through the subdomains - if mat.solubility_law == "henry": - F_henry += 1 * test_function_henry * mesh.dx(mat_id) - elif mat.solubility_law == "sievert": - F_sievert += 1 * test_function_sievert * mesh.dx(mat_id) - - # solve the problems - f.solve(F_henry == 0, henry, []) - f.solve(F_sievert == 0, sievert, []) - - self.henry_marker = henry - self.sievert_marker = sievert - - -class ArheniusCoeff(f.UserExpression): - def __init__(self, materials, vm, T, pre_exp, E, **kwargs): - super().__init__(kwargs) - self._vm = vm - self._T = T - self._materials = materials - self._pre_exp = pre_exp - self._E = E - - def eval_cell(self, value, x, ufc_cell): - cell = f.Cell(self._vm.mesh(), ufc_cell.index) - subdomain_id = self._vm[cell] - material = self._materials.find_material_from_id(subdomain_id) - D_0 = getattr(material, self._pre_exp) - E_D = getattr(material, self._E) - value[0] = D_0 * f.exp(-E_D / k_B / self._T(x)) - - def value_shape(self): - return () - - -class ThermalProp(f.UserExpression): - def __init__(self, materials, vm, T, key, **kwargs): - super().__init__(kwargs) - self._T = T - self._vm = vm - self._materials = materials - self._key = key - - def eval_cell(self, value, x, ufc_cell): - cell = f.Cell(self._vm.mesh(), ufc_cell.index) - subdomain_id = self._vm[cell] - material = self._materials.find_material_from_id(subdomain_id) - attribute = getattr(material, self._key) - if callable(attribute): - value[0] = attribute(self._T(x)) - else: - value[0] = attribute - - def value_shape(self): - return () diff --git a/meshing/mesh.py b/meshing/mesh.py deleted file mode 100644 index 27efb1fff..000000000 --- a/meshing/mesh.py +++ /dev/null @@ -1,37 +0,0 @@ -import fenics as f - - -class Mesh: - """ - Mesh class - - Args: - mesh (fenics.Mesh, optional): the mesh. Defaults to None. - volume_markers (fenics.MeshFunction, optional): markers of the mesh cells. Defaults to None. - surface_markers (fenics.MeshFunction, optional): markers of the mesh facets. Defaults to None. - type (str, optional): "cartesian", "cylindrical" or "spherical". Defaults to "cartesian". - - Attributes: - mesh (fenics.Mesh): the mesh - volume_markers (fenics.MeshFunction): markers of the mesh cells - surface_markers (fenics.MeshFunction): markers of the mesh facets - dx (fenics.Measure): - ds (fenics.Measure): - """ - - def __init__( - self, mesh=None, volume_markers=None, surface_markers=None, type="cartesian" - ) -> None: - self.mesh = mesh - self.volume_markers = volume_markers - self.surface_markers = surface_markers - self.type = type - - self.dx = None - self.ds = None - - def define_measures(self): - """Creates the fenics.Measure objects for self.dx and self.ds""" - - self.ds = f.Measure("ds", domain=self.mesh, subdomain_data=self.surface_markers) - self.dx = f.Measure("dx", domain=self.mesh, subdomain_data=self.volume_markers) diff --git a/meshing/mesh_1d.py b/meshing/mesh_1d.py deleted file mode 100644 index f845b8972..000000000 --- a/meshing/mesh_1d.py +++ /dev/null @@ -1,79 +0,0 @@ -import fenics as f -from festim import Mesh - - -class Mesh1D(Mesh): - """ - 1D Mesh - - Attributes: - size (float): the size of the 1D mesh - start (float): the starting point of the 1D mesh - - """ - - def __init__(self, **kwargs) -> None: - super().__init__(**kwargs) - self.size = None - self.start = 0 - - def define_markers(self, materials): - """Iterates through the mesh and mark them - based on their position in the domain - - Arguments: - materials {festim.Materials} -- contains the materials - """ - self.volume_markers = self.define_volume_markers(materials) - - self.surface_markers = self.define_surface_markers() - - def define_surface_markers(self): - """Creates the surface markers - - Returns: - fenics.MeshFunction: the meshfunction containing the surface - markers - """ - surface_markers = f.MeshFunction( - "size_t", self.mesh, self.mesh.topology().dim() - 1, 0 - ) - surface_markers.set_all(0) - i = 0 - for facet in f.facets(self.mesh): - i += 1 - x0 = facet.midpoint() - surface_markers[facet] = 0 - if f.near(x0.x(), self.start): - surface_markers[facet] = 1 - if f.near(x0.x(), self.size): - surface_markers[facet] = 2 - return surface_markers - - def define_volume_markers(self, materials): - """Creates the volume markers - - Args: - materials (festim.Materials): the materials - - Returns: - fenics.MeshFunction: the meshfunction containing the volume - markers - """ - volume_markers = f.MeshFunction( - "size_t", self.mesh, self.mesh.topology().dim(), 0 - ) - # iterate through the cells of the mesh and mark them - for cell in f.cells(self.mesh): - x = cell.midpoint().x() - subdomain_id = materials.find_subdomain_from_x_coordinate(x) - volume_markers[cell] = subdomain_id - - return volume_markers - - def define_measures(self, materials): - """Creates the fenics.Measure objects for self.dx and self.ds""" - if materials[0].borders is not None: - materials.check_borders(self.size) - self.define_markers(materials) - super().define_measures() diff --git a/meshing/mesh_from_vertices.py b/meshing/mesh_from_vertices.py deleted file mode 100644 index 7af680fe8..000000000 --- a/meshing/mesh_from_vertices.py +++ /dev/null @@ -1,44 +0,0 @@ -import fenics as f -import numpy as np -from festim import Mesh1D - - -class MeshFromVertices(Mesh1D): - """ - Description of MeshFromVertices - - Args: - vertices (list): the mesh vertices - - Attributes: - vertices (list): the mesh vertices - size (type): the size of the 1D mesh - """ - - def __init__(self, vertices, **kwargs) -> None: - super().__init__(**kwargs) - self.vertices = vertices - self.size = max(vertices) - self.start = min(vertices) - self.generate_mesh_from_vertices() - - def generate_mesh_from_vertices(self): - """Generates a 1D mesh""" - mesh = f.Mesh() - # the following is needed to avoid breaking in parrallel - # see issue 497 - if f.MPI.comm_world.rank == 0: - vertices = sorted(np.unique(self.vertices)) - nb_points = len(vertices) - nb_cells = nb_points - 1 - editor = f.MeshEditor() - editor.open(mesh, "interval", 1, 1) # top. and geom. dimension are both 1 - editor.init_vertices(nb_points) # number of vertices - editor.init_cells(nb_cells) # number of cells - for i in range(0, nb_points): - editor.add_vertex(i, np.array([vertices[i]])) - for j in range(0, nb_cells): - editor.add_cell(j, np.array([j, j + 1])) - editor.close() - f.MeshPartitioning.build_distributed_mesh(mesh) - self.mesh = mesh diff --git a/nonlinear_problem.py b/nonlinear_problem.py deleted file mode 100644 index f06bbd47e..000000000 --- a/nonlinear_problem.py +++ /dev/null @@ -1,31 +0,0 @@ -import fenics as f - - -class Problem(f.NonlinearProblem): - """ - Class to set up a nonlinear variational problem (F(u, v)=0) to solve - by the Newton method based on the form of the variational problem, the Jacobian - form of the variational problem, and the boundary conditions - - Args: - J (ufl.Form): the Jacobian form of the variational problem - F (ufl.Form): the form of the variational problem - bcs (list): list of fenics.DirichletBC - """ - - def __init__(self, J, F, bcs): - self.jacobian_form = J - self.residual_form = F - self.bcs = bcs - self.assembler = f.SystemAssembler( - self.jacobian_form, self.residual_form, self.bcs - ) - f.NonlinearProblem.__init__(self) - - def F(self, b, x): - """Assembles the RHS in Ax=b and applies the boundary conditions""" - self.assembler.assemble(b, x) - - def J(self, A, x): - """Assembles the LHS in Ax=b and applies the boundary conditions""" - self.assembler.assemble(A) diff --git a/settings.py b/settings.py deleted file mode 100644 index 533fba778..000000000 --- a/settings.py +++ /dev/null @@ -1,76 +0,0 @@ -class Settings: - """ - Args: - absolute_tolerance (float): the absolute tolerance of the newton - solver - relative_tolerance (float): the relative tolerance of the newton - solver - maximum_iterations (int, optional): maximum iterations allowed for - the solver to converge. Defaults to 30. - transient (bool, optional): If set to True, the simulation will be - transient. Defaults to True. - final_time (float, optional): final time of the simulation. - Defaults to None. - chemical_pot (bool, optional): if True, conservation of chemical - potential will be assumed. Defaults to False. - soret (bool, optional): if True, Soret effect will be assumed. - Defaults to False. - traps_element_type (str, optional): Finite element used for traps. - If traps densities are discontinuous (eg. different materials) - "DG" is recommended. Defaults to "CG". - update_jacobian (bool, optional): If set to False, the Jacobian of - the formulation will be computed only once at the beggining. - Else it will be computed at each time step. Defaults to True. - linear_solver (str, optional): linear solver method for the newton solver, - options can be viewed by print(list_linear_solver_methods()). - More information can be found at: https://fenicsproject.org/pub/tutorial/html/._ftut1017.html. - Defaults to None, for the newton solver this is: "umfpack". - preconditioner (str, optional): preconditioning method for the newton solver, - options can be viewed by print(list_krylov_solver_preconditioners()). - Defaults to "default". - - Attributes: - transient (bool): transient or steady state sim - final_time (float): final time of the simulation. - chemical_pot (bool): conservation of chemical potential - soret (bool): soret effect - absolute_tolerance (float): the absolute tolerance of the newton - solver - relative_tolerance (float): the relative tolerance of the newton - solver - maximum_iterations (int): maximum iterations allowed for - the solver to converge - traps_element_type (str): Finite element used for traps. - update_jacobian (bool): - linear_solver (str): linear solver method for the newton solver - precondtitioner (str): preconditioning method for the newton solver - """ - - def __init__( - self, - absolute_tolerance, - relative_tolerance, - maximum_iterations=30, - transient=True, - final_time=None, - chemical_pot=False, - soret=False, - traps_element_type="CG", - update_jacobian=True, - linear_solver=None, - preconditioner="default", - ): - # TODO maybe transient and final_time are redundant - self.transient = transient - self.final_time = final_time - self.chemical_pot = chemical_pot - self.soret = soret - - self.absolute_tolerance = absolute_tolerance - self.relative_tolerance = relative_tolerance - self.maximum_iterations = maximum_iterations - - self.traps_element_type = traps_element_type - self.update_jacobian = update_jacobian - self.linear_solver = linear_solver - self.preconditioner = preconditioner diff --git a/sources/radioactive_decay.py b/sources/radioactive_decay.py deleted file mode 100644 index 6150a3dcd..000000000 --- a/sources/radioactive_decay.py +++ /dev/null @@ -1,35 +0,0 @@ -from festim import Source - - -class RadioactiveDecay(Source): - """ - A radioactive decay source. - dc/dt = ... - lambda * c - where lambda is the decay constant in 1/s. - - Args: - decay_constant (float): The decay constant of the source in 1/s. - volume (int): The volume of the source. - field (str, optional): The field to which the source is - applied. If "all" the decay will be applied to all - concentrations. Defaults to "all". - """ - - def __init__(self, decay_constant, volume, field="all") -> None: - self.decay_constant = decay_constant - super().__init__(value=None, volume=volume, field=field) - - @property - def decay_constant(self): - return self._decay_constant - - @decay_constant.setter - def decay_constant(self, value): - if not isinstance(value, (float, int)): - raise TypeError("decay_constant must be a float or an int") - if value <= 0: - raise ValueError("decay_constant must be positive") - self._decay_constant = value - - def form(self, concentration): - return -self.decay_constant * concentration diff --git a/stepsize.py b/stepsize.py deleted file mode 100644 index f60106904..000000000 --- a/stepsize.py +++ /dev/null @@ -1,139 +0,0 @@ -import fenics as f -import numpy as np -import warnings - - -class Stepsize: - """ - Description of Stepsize - - Args: - initial_value (float, optional): initial stepsize. Defaults to 0.0. - stepsize_change_ratio (float, optional): stepsize change ratio. - Defaults to None. - t_stop (float, optional): time at which the adaptive stepsize - stops. Defaults to None. - stepsize_stop_max (float, optional): Maximum stepsize after - t_stop. Defaults to None. - max_stepsize (float or callable, optional): Maximum stepsize. - Can be a function of festim.t. Defaults to None. - dt_min (float, optional): Minimum stepsize below which an error is - raised. Defaults to None. - milestones (list, optional): list of times by which the simulation must - pass. Defaults to None. - - Attributes: - adaptive_stepsize (dict): contains the parameters for adaptive stepsize - value (fenics.Constant): value of dt - milestones (list): list of times by which the simulation must - pass. - - Example:: - - my_stepsize = Stepsize( - initial_value=0.5, - stepsize_change_ratio=1.1, - max_stepsize=lambda t: None if t < 1 else 2, - dt_min=1e-05 - ) - """ - - def __init__( - self, - initial_value=0.0, - stepsize_change_ratio=None, - t_stop=None, - stepsize_stop_max=None, - max_stepsize=None, - dt_min=None, - milestones=None, - ) -> None: - self.adaptive_stepsize = None - if stepsize_change_ratio is not None: - if t_stop or stepsize_stop_max: - warnings.warn( - "stepsize_stop_max and t_stop attributes will be deprecated in a future release, please use max_stepsize instead", - DeprecationWarning, - ) - max_stepsize = lambda t: stepsize_stop_max if t >= t_stop else None - self.adaptive_stepsize = { - "stepsize_change_ratio": stepsize_change_ratio, - "max_stepsize": max_stepsize, - "dt_min": dt_min, - } - self.initial_value = initial_value - self.value = None - self.milestones = milestones - self.initialise_value() - - @property - def milestones(self): - return self._milestones - - @milestones.setter - def milestones(self, value): - if value: - self._milestones = sorted(value) - else: - self._milestones = value - - def initialise_value(self): - """Creates a fenics.Constant object initialised with self.initial_value - and stores it in self.value""" - self.value = f.Constant(self.initial_value, name="dt") - - def adapt(self, t, nb_it, converged): - """Changes the stepsize based on convergence. - - Args: - t (float): current time. - nb_it (int): number of iterations the solver required to converge. - converged (bool): True if the solver converged, else False. - """ - if self.adaptive_stepsize: - change_ratio = self.adaptive_stepsize["stepsize_change_ratio"] - dt_min = self.adaptive_stepsize["dt_min"] - max_stepsize = self.adaptive_stepsize["max_stepsize"] - - if not converged: - if dt_min is None: - raise ValueError("Solver diverged but dt_min is not set.") - - self.value.assign(float(self.value) / change_ratio) - if float(self.value) < dt_min: - raise ValueError("stepsize reached minimal value") - if nb_it < 5: - self.value.assign(float(self.value) * change_ratio) - else: - self.value.assign(float(self.value) / change_ratio) - - if callable(max_stepsize): - max_stepsize = max_stepsize(t) - if max_stepsize is not None: - if float(self.value) > max_stepsize: - self.value.assign(max_stepsize) - - # adapt for next milestone - next_milestone = self.next_milestone(t) - if next_milestone is not None: - if t + float(self.value) > next_milestone and not np.isclose( - t, next_milestone, atol=0 - ): - self.value.assign((next_milestone - t)) - - def next_milestone(self, current_time: float): - """Returns the next milestone that the simulation must pass. - Returns None if there are no more milestones. - - Args: - current_time (float): current time. - - Returns: - float: next milestone. - """ - if self.milestones is None: - return None - for milestone in self.milestones: - if current_time < milestone: - return milestone - return None diff --git a/temperature/temperature.py b/temperature/temperature.py deleted file mode 100644 index be695c772..000000000 --- a/temperature/temperature.py +++ /dev/null @@ -1,56 +0,0 @@ -import sympy as sp -import fenics as f - - -class Temperature: - """ - Class for Temperature in FESTIM - - Args: - value (sp.Add, int, float, optional): The value of the temperature. - Defaults to None. - - Attributes: - T (fenics.Function): the function attributed with temperature - T_n (fenics.Function): the previous function - value (sp.Add, int, float): the expression of temperature - expression (fenics.Expression): the expression of temperature as a - fenics object - - Usage: - >>> import festim as F - >>> my_model = F.Simulation(...) - >>> my_model.T = F.Temperature(300 + 10 * F.x + F.t) - """ - - def __init__(self, value=None) -> None: - self.T = None - self.T_n = None - self.value = value - self.expression = None - - def create_functions(self, mesh): - """Creates functions self.T, self.T_n - - Args: - mesh (festim.Mesh): the mesh - """ - V = f.FunctionSpace(mesh.mesh, "CG", 1) - self.T = f.Function(V, name="T") - self.T_n = f.Function(V, name="T_n") - self.expression = f.Expression(sp.printing.ccode(self.value), t=0, degree=2) - self.T.assign(f.interpolate(self.expression, V)) - self.T_n.assign(self.T) - - def update(self, t): - """Updates T_n, expression, and T with respect to time - - Args: - t (float): the time - """ - self.T_n.assign(self.T) - self.expression.t = t - self.T.assign(f.interpolate(self.expression, self.T.function_space())) - - def is_steady_state(self): - return "t" not in sp.printing.ccode(self.value) diff --git a/temperature/temperature_from_xdmf.py b/temperature/temperature_from_xdmf.py deleted file mode 100644 index 58bd77338..000000000 --- a/temperature/temperature_from_xdmf.py +++ /dev/null @@ -1,52 +0,0 @@ -from festim.temperature.temperature import Temperature -from festim.helpers import extract_xdmf_labels -import fenics as f - - -class TemperatureFromXDMF(Temperature): - """ - Temperature read from an XDMF file - - Args: - filename (str): The temperature file. Must end in ".xdmf" - label (str): How the checkpoints have been labelled - - Attributes: - filename (str): name of the temperature file - label (str): How the checkpoints have been labelled - """ - - def __init__(self, filename, label) -> None: - super().__init__() - - self.filename = filename - self.label = label - - # check labels match - if self.label not in extract_xdmf_labels(self.filename): - raise ValueError( - "Coudln't find label: {} in {}".format(self.label, self.filename) - ) - - def create_functions(self, mesh): - """Creates functions self.T, self.T_n - Args: - mesh (festim.Mesh): the mesh - """ - V = f.FunctionSpace(mesh.mesh, "CG", 1) - self.T = f.Function(V, name="T") - - f.XDMFFile(self.filename).read_checkpoint(self.T, self.label, -1) - - self.T_n = f.Function(V, name="T_n") - self.T_n.assign(self.T) - - def update(self, t): - """Allows for the use of this class in transient h transport cases, - refer to issue #499 - """ - pass - - def is_steady_state(self): - # TemperatureFromXDMF is always steady state - return True diff --git a/temperature/temperature_solver.py b/temperature/temperature_solver.py deleted file mode 100644 index c35a635a0..000000000 --- a/temperature/temperature_solver.py +++ /dev/null @@ -1,275 +0,0 @@ -import festim -import fenics as f -import sympy as sp - - -class HeatTransferProblem(festim.Temperature): - """ - Args: - transient (bool, optional): If True, a transient simulation will - be run. Defaults to True. - initial_condition (int, float, sp.Expr, festim.InitialCondition, optional): The initial condition. - Only needed if transient is True. - absolute_tolerance (float, optional): the absolute tolerance of the newton - solver. Defaults to 1e-03 - relative_tolerance (float, optional): the relative tolerance of the newton - solver. Defaults to 1e-10 - maximum_iterations (int, optional): maximum iterations allowed for - the solver to converge. Defaults to 30. - linear_solver (str, optional): linear solver method for the newton solver, - options can be viewed with print(list_linear_solver_methods()). - If None, the default fenics linear solver will be used ("umfpack"). - More information can be found at: https://fenicsproject.org/pub/tutorial/html/._ftut1017.html. - Defaults to None. - preconditioner (str, optional): preconditioning method for the newton solver, - options can be veiwed by print(list_krylov_solver_preconditioners()). - Defaults to "default". - - Attributes: - F (fenics.Form): the variational form of the heat transfer problem - v_T (fenics.TestFunction): the test function - newton_solver (fenics.NewtonSolver): Newton solver for solving the nonlinear problem - initial_condition (festim.InitialCondition): the initial condition - sub_expressions (list): contains time dependent fenics.Expression to - be updated - sources (list): contains festim.Source objects for volumetric heat - sources - boundary_conditions (list): contains festim.BoundaryConditions - """ - - def __init__( - self, - transient=True, - initial_condition=None, - absolute_tolerance=1e-3, - relative_tolerance=1e-10, - maximum_iterations=30, - linear_solver=None, - preconditioner="default", - ) -> None: - super().__init__() - self.transient = transient - self.initial_condition = initial_condition - self.absolute_tolerance = absolute_tolerance - self.relative_tolerance = relative_tolerance - self.maximum_iterations = maximum_iterations - self.linear_solver = linear_solver - self.preconditioner = preconditioner - - self.F = 0 - self.v_T = None - self.sources = [] - self.boundary_conditions = [] - self.sub_expressions = [] - self.newton_solver = None - - @property - def newton_solver(self): - return self._newton_solver - - @newton_solver.setter - def newton_solver(self, value): - if value is None: - self._newton_solver = value - elif isinstance(value, f.NewtonSolver): - if self._newton_solver: - print("Settings for the Newton solver will be overwritten") - self._newton_solver = value - else: - raise TypeError("accepted type for newton_solver is fenics.NewtonSolver") - - @property - def initial_condition(self): - return self._initial_condition - - @initial_condition.setter - def initial_condition(self, value): - if isinstance(value, (int, float, sp.Expr)): - self._initial_condition = festim.InitialCondition(field="T", value=value) - else: - self._initial_condition = value - - # TODO rename initialise? - def create_functions(self, materials, mesh, dt=None): - """Creates functions self.T, self.T_n and test function self.v_T. - Solves the steady-state heat transfer problem if self.transient is - False. - - Args: - materials (festim.Materials): the materials. - mesh (festim.Mesh): the mesh - dt (festim.Stepsize, optional): the stepsize. Only needed if - self.transient is True. Defaults to None. - """ - # Define variational problem for heat transfers - V = f.FunctionSpace(mesh.mesh, "CG", 1) - self.T = f.Function(V, name="T") - self.T_n = f.Function(V, name="T_n") - self.v_T = f.TestFunction(V) - if self.transient and self.initial_condition is None: - raise AttributeError( - "Initial condition is required for transient heat transfer simulations" - ) - if self.transient and self.initial_condition: - if isinstance(self.initial_condition.value, str): - if self.initial_condition.value.endswith(".xdmf"): - with f.XDMFFile(self.initial_condition.value) as file: - file.read_checkpoint( - self.T_n, - self.initial_condition.label, - self.initial_condition.time_step, - ) - else: - ccode_T_ini = sp.printing.ccode(self.initial_condition.value) - self.initial_condition.value = f.Expression(ccode_T_ini, degree=2, t=0) - self.T_n.assign(f.interpolate(self.initial_condition.value, V)) - - self.define_variational_problem(materials, mesh, dt) - self.create_dirichlet_bcs(mesh.surface_markers) - - if not self.newton_solver: - self.define_newton_solver() - - if not self.transient: - print("Solving stationary heat equation") - dT = f.TrialFunction(self.T.function_space()) - JT = f.derivative(self.F, self.T, dT) - problem = festim.Problem(JT, self.F, self.dirichlet_bcs) - - f.begin( - "Solving nonlinear variational problem." - ) # Add message to fenics logs - self.newton_solver.solve(problem, self.T.vector()) - f.end() - - self.T_n.assign(self.T) - - def define_variational_problem(self, materials, mesh, dt=None): - """Create a variational form for heat transfer problem - - Args: - materials (festim.Materials): the materials. - mesh (festim.Mesh): the mesh. - dt (festim.Stepsize, optional): the stepsize. Only needed if - self.transient is True. Defaults to None. - """ - - print("Defining variational problem heat transfers") - T, T_n = self.T, self.T_n - v_T = self.v_T - - self.F = 0 - for mat in materials: - thermal_cond = mat.thermal_cond - if callable(thermal_cond): # if thermal_cond is a function - thermal_cond = thermal_cond(T) - - subdomains = mat.id # list of subdomains with this material - if type(subdomains) is not list: - subdomains = [subdomains] # make sure subdomains is a list - if self.transient: - cp = mat.heat_capacity - rho = mat.rho - if callable(cp): # if cp or rho are functions, apply T - cp = cp(T) - if callable(rho): - rho = rho(T) - # Transien term - for vol in subdomains: - self.F += rho * cp * (T - T_n) / dt.value * v_T * mesh.dx(vol) - # Diffusion term - for vol in subdomains: - if mesh.type == "cartesian": - self.F += f.dot(thermal_cond * f.grad(T), f.grad(v_T)) * mesh.dx( - vol - ) - elif mesh.type == "cylindrical": - r = f.SpatialCoordinate(mesh.mesh)[0] - self.F += ( - r - * f.dot(thermal_cond * f.grad(T), f.grad(v_T / r)) - * mesh.dx(vol) - ) - elif mesh.type == "spherical": - r = f.SpatialCoordinate(mesh.mesh)[0] - self.F += ( - thermal_cond - * r - * r - * f.dot(f.grad(T), f.grad(v_T / r / r)) - * mesh.dx(vol) - ) - # source term - for source in self.sources: - self.sub_expressions.append(source.value) - if type(source.volume) is list: - volumes = source.volume - else: - volumes = [source.volume] - for volume in volumes: - self.F += -source.value * v_T * mesh.dx(volume) - - # Boundary conditions - for bc in self.boundary_conditions: - if isinstance(bc, festim.FluxBC): - bc.create_form(self.T, solute=None) - - # TODO: maybe that's not necessary - self.sub_expressions += bc.sub_expressions - - for surf in bc.surfaces: - self.F += -bc.form * self.v_T * mesh.ds(surf) - - def define_newton_solver(self): - """Creates the Newton solver and sets its parameters""" - self.newton_solver = f.NewtonSolver(f.MPI.comm_world) - self.newton_solver.parameters["error_on_nonconvergence"] = True - self.newton_solver.parameters["absolute_tolerance"] = self.absolute_tolerance - self.newton_solver.parameters["relative_tolerance"] = self.relative_tolerance - self.newton_solver.parameters["maximum_iterations"] = self.maximum_iterations - self.newton_solver.parameters["linear_solver"] = self.linear_solver - self.newton_solver.parameters["preconditioner"] = self.preconditioner - - def create_dirichlet_bcs(self, surface_markers): - """Creates a list of fenics.DirichletBC and add time dependent - expressions to .sub_expressions - - Args: - surface_markers (fenics.MeshFunction): contains the mesh facet - markers - """ - V = self.T.function_space() - self.dirichlet_bcs = [] - for bc in self.boundary_conditions: - if isinstance(bc, festim.DirichletBC) and bc.field == "T": - bc.create_expression(self.T) - for surf in bc.surfaces: - bci = f.DirichletBC(V, bc.expression, surface_markers, surf) - self.dirichlet_bcs.append(bci) - self.sub_expressions += bc.sub_expressions - self.sub_expressions.append(bc.expression) - - def update(self, t): - """Updates T_n, and T with respect to time by solving the heat transfer - problem - - Args: - t (float): the time - """ - if self.transient: - festim.update_expressions(self.sub_expressions, t) - # Solve heat transfers - dT = f.TrialFunction(self.T.function_space()) - JT = f.derivative(self.F, self.T, dT) # Define the Jacobian - problem = festim.Problem(JT, self.F, self.dirichlet_bcs) - - f.begin( - "Solving nonlinear variational problem." - ) # Add message to fenics logs - self.newton_solver.solve(problem, self.T.vector()) - f.end() - - self.T_n.assign(self.T) - - def is_steady_state(self): - return not self.transient