iterate method

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[:]