Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Settings and Stepsize classes #616

Merged
merged 21 commits into from
Oct 20, 2023
Merged
Show file tree
Hide file tree
Changes from 11 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions festim/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -28,10 +28,14 @@

from .hydrogen_transport_problem import HydrogenTransportProblem

from .settings import Settings

from .species import Species, Trap

from .subdomain.surface_subdomain import SurfaceSubdomain1D
from .subdomain.volume_subdomain import VolumeSubdomain1D

from .stepsize import Stepsize

from .exports.vtx import VTXExport
from .exports.xdmf import XDMFExport
30 changes: 15 additions & 15 deletions festim/hydrogen_transport_problem.py
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,7 @@ def __init__(
temperature=None,
sources=[],
boundary_conditions=[],
solver_parameters=None,
settings=None,
exports=[],
) -> None:
self.mesh = mesh
Expand All @@ -86,7 +86,7 @@ def __init__(
self.temperature = temperature
self.sources = sources
self.boundary_conditions = boundary_conditions
self.solver_parameters = solver_parameters
self.settings = settings
self.exports = exports

self.dx = None
Expand Down Expand Up @@ -218,10 +218,9 @@ def create_formulation(self):
if len(self.species) > 1:
raise NotImplementedError("Multiple species not implemented yet")

# TODO expose dt as parameter of the model
dt = fem.Constant(self.mesh.mesh, 1 / 20)

self.dt = dt # TODO remove this
self.dt = F.as_fenics_constant(
self.settings.stepsize.initial_value, self.mesh.mesh
)
jhdark marked this conversation as resolved.
Show resolved Hide resolved

self.formulation = 0

Expand All @@ -236,7 +235,7 @@ def create_formulation(self):
)

self.formulation += dot(D * grad(u), grad(v)) * self.dx(vol.id)
self.formulation += ((u - u_n) / dt) * v * self.dx(vol.id)
self.formulation += ((u - u_n) / self.dt) * v * self.dx(vol.id)

# add sources
# TODO implement this
Expand All @@ -258,15 +257,14 @@ def create_solver(self):
self.species[0].solution,
bcs=self.bc_forms,
)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
self.solver = solver
self.solver = NewtonSolver(MPI.COMM_WORLD, problem)
self.solver.atol = self.settings.absolute_tolerance
self.solver.rtol = self.settings.relative_tolerance
self.solver.max_it = self.settings.max_iterations
jhdark marked this conversation as resolved.
Show resolved Hide resolved

def run(self, final_time: float):
def run(self):
"""Runs the model for a given time

Args:
final_time (float): the final time of the simulation

Returns:
list of float: the times of the simulation
list of float: the fluxes of the simulation
Expand All @@ -279,9 +277,11 @@ def run(self, final_time: float):
)
cm = self.species[0].solution
progress = tqdm.autonotebook.tqdm(
desc="Solving H transport problem", total=final_time
desc="Solving H transport problem",
total=self.settings.final_time,
unit_scale=True,
jhdark marked this conversation as resolved.
Show resolved Hide resolved
)
while self.t.value < final_time:
while self.t.value < self.settings.final_time:
progress.update(self.dt.value)
self.t.value += self.dt.value

Expand Down
33 changes: 33 additions & 0 deletions festim/settings.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,33 @@
class Settings:
"""Settings for a festim simulation.

Args:
absolute_tolerance (float): Absolute tolerance for the solver.
relative_tolerance (float): Relative tolerance for the solver.
max_iterations (int, optional): Maximum number of iterations for the solver.
final_time (float, optional): Final time for a transient simulation.
stepsize (festim.Stepsize, optional): stepsize for a transient
simulation.

Attributes:
aboslute_tolerance (float): Absolute tolerance for the solver.
jhdark marked this conversation as resolved.
Show resolved Hide resolved
relative_tolerance (float): Relative tolerance for the solver.
max_iterations (int, optional): Maximum number of iterations for the solver.
jhdark marked this conversation as resolved.
Show resolved Hide resolved
final_time (float, optional): Final time for a transient simulation.
stepsize (festim.Stepsize, optional): stepsize for a transient
jhdark marked this conversation as resolved.
Show resolved Hide resolved
jhdark marked this conversation as resolved.
Show resolved Hide resolved
simulation.
"""

def __init__(
self,
absolute_tolerance,
relative_tolerance,
jhdark marked this conversation as resolved.
Show resolved Hide resolved
max_iterations=30,
final_time=None,
stepsize=None,
) -> None:
self.absolute_tolerance = absolute_tolerance
self.relative_tolerance = relative_tolerance
self.max_iterations = max_iterations
self.final_time = final_time
self.stepsize = stepsize
16 changes: 16 additions & 0 deletions festim/stepsize.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,16 @@
class Stepsize:
"""
A class for evaulating the stepsize of transient simulations.
jhdark marked this conversation as resolved.
Show resolved Hide resolved

Args:
initial_value (float, int): initial stepsize.

Attributes:
initial_value (float, int): initial stepsize.
"""

def __init__(
self,
initial_value,
) -> None:
self.initial_value = initial_value
7 changes: 6 additions & 1 deletion test/test_dirichlet_bc.py
Original file line number Diff line number Diff line change
Expand Up @@ -281,13 +281,18 @@ def test_integration_with_HTransportProblem(value):

my_model.temperature = fem.Constant(my_model.mesh.mesh, 550.0)

my_model.settings = F.Settings(
absolute_tolerance=1, relative_tolerance=0.1, final_time=2
)
my_model.settings.stepsize = F.Stepsize(initial_value=1)

# RUN

my_model.initialise()

assert my_bc.value_fenics is not None

my_model.run(final_time=2)
my_model.run()

# TEST

Expand Down
34 changes: 20 additions & 14 deletions test/test_permeation_problem.py
Original file line number Diff line number Diff line change
Expand Up @@ -35,13 +35,18 @@ def test_permeation_problem(mesh_size=1001):
]
my_model.exports = [F.XDMFExport("mobile_concentration.xdmf", field=mobile_H)]

my_model.settings = F.Settings(
absolute_tolerance=1e10,
relative_tolerance=1e-10,
max_iterations=30,
final_time=50,
)

my_model.settings.stepsize = F.Stepsize(initial_value=1 / 20)

my_model.initialise()

my_model.solver.convergence_criterion = "incremental"
my_model.solver.rtol = 1e-10
my_model.solver.atol = 1e10

my_model.solver.report = True
ksp = my_model.solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
Expand All @@ -50,9 +55,7 @@ def test_permeation_problem(mesh_size=1001):
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()

final_time = 50

times, flux_values = my_model.run(final_time=final_time)
times, flux_values = my_model.run()

# -------------------------- analytical solution -------------------------------------

Expand Down Expand Up @@ -126,13 +129,18 @@ def test_permeation_problem_multi_volume():
]
my_model.exports = [F.VTXExport("test.bp", field=mobile_H)]

my_model.settings = F.Settings(
absolute_tolerance=1e10,
relative_tolerance=1e-10,
max_iterations=30,
final_time=50,
)

my_model.settings.stepsize = F.Stepsize(initial_value=1 / 20)

my_model.initialise()

my_model.solver.convergence_criterion = "incremental"
my_model.solver.rtol = 1e-10
my_model.solver.atol = 1e10

my_model.solver.report = True
ksp = my_model.solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
Expand All @@ -141,9 +149,7 @@ def test_permeation_problem_multi_volume():
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()

final_time = 50

times, flux_values = my_model.run(final_time=final_time)
times, flux_values = my_model.run()

# -------------------------- analytical solution -------------------------------------
D = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature)
Expand Down
7 changes: 6 additions & 1 deletion test/test_sievertsbc.py
Original file line number Diff line number Diff line change
Expand Up @@ -93,10 +93,15 @@ def test_integration_with_HTransportProblem(pressure):

my_model.temperature = fem.Constant(my_model.mesh.mesh, 550.0)

my_model.settings = F.Settings(
absolute_tolerance=1, relative_tolerance=0.1, final_time=2
)
my_model.settings.stepsize = F.Stepsize(initial_value=1)

# RUN

my_model.initialise()

assert my_bc.value_fenics is not None

my_model.run(final_time=2)
my_model.run()
2 changes: 2 additions & 0 deletions test/test_vtx.py
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,8 @@ def test_vtx_integration_with_h_transport_problem(tmpdir):
filename = str(tmpdir.join("my_export.bp"))
my_export = F.VTXExport(filename, field=my_model.species[0])
my_model.exports = [my_export]
my_model.settings = F.Settings(absolute_tolerance=1, relative_tolerance=0.1)
my_model.settings.stepsize = F.Stepsize(initial_value=1)

my_model.initialise()

Expand Down
6 changes: 5 additions & 1 deletion test/test_xdmf.py
Original file line number Diff line number Diff line change
Expand Up @@ -47,8 +47,12 @@ def test_integration_with_HTransportProblem(tmp_path):
filename = os.path.join(tmp_path, "test.xdmf")
my_model.exports = [F.XDMFExport(filename=filename, field=my_model.species)]

my_model.settings = F.Settings(
absolute_tolerance=1, relative_tolerance=0.1, final_time=1
)
my_model.settings.stepsize = F.Stepsize(initial_value=0.5)
my_model.initialise()
my_model.run(1)
my_model.run()

# checks that filename exists
assert os.path.exists(filename)
Expand Down
Loading