From 6229d73d9e259aae4a3233fdede49117ed2cdd81 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 31 Oct 2023 08:41:40 -0400 Subject: [PATCH] iterate method --- festim/hydrogen_transport_problem.py | 98 +++++++++++++++------------- 1 file changed, 51 insertions(+), 47 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index ff171bca3..61b8de418 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -357,70 +357,74 @@ def create_solver(self): self.solver.max_it = self.settings.max_iterations def run(self): - """Runs the model for a given time + """Runs the model Returns: list of float: the times of the simulation list of float: the fluxes of the simulation """ - times, flux_values = [], [] - flux_values_1, flux_values_2 = [], [] + self.times, self.flux_values = [], [] + self.flux_values_1, self.flux_values_2 = [], [] - progress = tqdm.autonotebook.tqdm( + self.progress = tqdm.autonotebook.tqdm( desc="Solving H transport problem", total=self.settings.final_time, unit_scale=True, ) while self.t.value < self.settings.final_time: - progress.update(self.dt.value) - self.t.value += self.dt.value + self.iterate() - # update boundary conditions - for bc in self.boundary_conditions: - bc.update(float(self.t)) + if self.multispecies: + self.flux_values = [self.flux_values_1, self.flux_values_2] - self.solver.solve(self.u) + return self.times, self.flux_values - # post processing - # TODO remove this - if not self.multispecies: - D_D = self.subdomains[0].material.get_diffusion_coefficient( - self.mesh.mesh, self.temperature, self.species[0] - ) - cm = self.u - self.species[0].post_processing_solution = self.u + def iterate(self): + """Iterates the model for a given time step""" + self.progress.update(self.dt.value) + self.t.value += self.dt.value - surface_flux = form(D_D * dot(grad(cm), self.mesh.n) * self.ds(2)) - flux = assemble_scalar(surface_flux) - flux_values.append(flux) - times.append(float(self.t)) - else: - for idx, spe in enumerate(self.species): - spe.post_processing_solution = self.u.sub(idx) + # update boundary conditions + for bc in self.boundary_conditions: + bc.update(float(self.t)) - cm_1, cm_2 = self.u.split() - D_1 = self.subdomains[0].material.get_diffusion_coefficient( - self.mesh.mesh, self.temperature, self.species[0] - ) - D_2 = self.subdomains[0].material.get_diffusion_coefficient( - self.mesh.mesh, self.temperature, self.species[1] - ) - surface_flux_1 = form(D_1 * dot(grad(cm_1), self.mesh.n) * self.ds(2)) - surface_flux_2 = form(D_2 * dot(grad(cm_2), self.mesh.n) * self.ds(2)) - flux_1 = assemble_scalar(surface_flux_1) - flux_2 = assemble_scalar(surface_flux_2) - flux_values_1.append(flux_1) - flux_values_2.append(flux_2) - times.append(float(self.t)) + self.solver.solve(self.u) - for export in self.exports: - if isinstance(export, (F.VTXExport, F.XDMFExport)): - export.write(float(self.t)) + # post processing + # TODO remove this + if not self.multispecies: + D_D = self.subdomains[0].material.get_diffusion_coefficient( + self.mesh.mesh, self.temperature, self.species[0] + ) + cm = self.u + self.species[0].post_processing_solution = self.u - # update previous solution - self.u_n.x.array[:] = self.u.x.array[:] + surface_flux = form(D_D * dot(grad(cm), self.mesh.n) * self.ds(2)) + flux = assemble_scalar(surface_flux) + self.flux_values.append(flux) + self.times.append(float(self.t)) + else: + for idx, spe in enumerate(self.species): + spe.post_processing_solution = self.u.sub(idx) - if self.multispecies: - flux_values = [flux_values_1, flux_values_2] + cm_1, cm_2 = self.u.split() + D_1 = self.subdomains[0].material.get_diffusion_coefficient( + self.mesh.mesh, self.temperature, self.species[0] + ) + D_2 = self.subdomains[0].material.get_diffusion_coefficient( + self.mesh.mesh, self.temperature, self.species[1] + ) + surface_flux_1 = form(D_1 * dot(grad(cm_1), self.mesh.n) * self.ds(2)) + surface_flux_2 = form(D_2 * dot(grad(cm_2), self.mesh.n) * self.ds(2)) + flux_1 = assemble_scalar(surface_flux_1) + flux_2 = assemble_scalar(surface_flux_2) + self.flux_values_1.append(flux_1) + self.flux_values_2.append(flux_2) + self.times.append(float(self.t)) + + for export in self.exports: + if isinstance(export, (F.VTXExport, F.XDMFExport)): + export.write(float(self.t)) - return times, flux_values + # update previous solution + self.u_n.x.array[:] = self.u.x.array[:]