Merge pull request #23 from KVSlab/store-u-and-d

Store mesh velocity and deformation during checkpointing.

codecov.yml

@@ -0,0 +1,2 @@
+ignore:
+  - "tests/*"

src/oasismove/NSfracStep.py

@@ -111,6 +111,9 @@
 q_1 = dict((ui, Function(VV[ui], name=ui + "_1")) for ui in sys_comp)
 q_2 = dict((ui, Function(V, name=ui + "_2")) for ui in u_components)
 
+# Create temporary dictionary for the mesh velocity and deformation
+w_ = d_ = None
+
 # Read in previous solution if restarting
 init_from_restart(**vars())
 
@@ -119,6 +122,7 @@
 u_1 = as_vector([q_1[ui] for ui in u_components])  # Velocity vector at t - dt
 u_2 = as_vector([q_2[ui] for ui in u_components])  # Velocity vector at t - 2*dt
 
+
 # Adams Bashforth projection of velocity at t - dt/2
 U_AB = 1.5 * u_1 - 0.5 * u_2
 

src/oasismove/common/io.py

@@ -9,6 +9,7 @@
 from xml.etree import ElementTree as ET
 
 from dolfin import MPI, XDMFFile, HDF5File, FunctionSpace, Function, interpolate, MeshFunction
+
 from oasismove.problems import info_red
 
 __all__ = ["create_initial_folders", "save_solution", "save_tstep_solution_xdmf",
@@ -62,8 +63,8 @@ def create_initial_folders(folder, restart_folder, sys_comp, tstep, info_red,
     return newfolder, tstepfiles
 
 
-def save_solution(tstep, t, q_, q_1, folder, newfolder, save_step, checkpoint, NS_parameters, tstepfiles, u_,
-                  u_components, output_timeseries_as_vector, mesh, AssignedVectorFunction, killtime, total_timer,
+def save_solution(tstep, t, q_, q_1, w_, d_, folder, newfolder, save_step, checkpoint, NS_parameters, tstepfiles,
+                  u_, u_components, output_timeseries_as_vector, mesh, AssignedVectorFunction, killtime, total_timer,
                   **NS_namespace):
     """Called at end of timestep. Check for kill and save solution if required."""
     NS_parameters.update(t=t, tstep=tstep)
@@ -78,7 +79,7 @@ def save_solution(tstep, t, q_, q_1, folder, newfolder, save_step, checkpoint, N
 
     killoasis = check_if_kill(folder, killtime, total_timer)
     if tstep % checkpoint == 0 or killoasis:
-        save_checkpoint_solution_xdmf(q_, q_1, newfolder, u_components, mesh, NS_parameters)
+        save_checkpoint_solution_xdmf(q_, q_1, w_, d_, newfolder, u_components, mesh, NS_parameters)
 
     return killoasis
 
@@ -116,7 +117,7 @@ def save_tstep_solution_xdmf(tstep, q_, u_, newfolder, tstepfiles, output_timese
             pickle.dump(NS_parameters, f)
 
 
-def save_checkpoint_solution_xdmf(q_, q_1, newfolder, u_components, mesh, NS_parameters):
+def save_checkpoint_solution_xdmf(q_, q_1, w_, d_, newfolder, u_components, mesh, NS_parameters):
     """
     Overwrite solution in Checkpoint folder
     """
@@ -142,6 +143,19 @@ def save_checkpoint_solution_xdmf(q_, q_1, newfolder, u_components, mesh, NS_par
                 f.write_checkpoint(q_1[ui], '/previous', append=True)
         MPI.barrier(MPI.comm_world)
 
+    MPI.barrier(MPI.comm_world)
+    # Store mesh velocity and deformation solution
+    if w_ is not None:
+        for ui in w_:
+            checkpoint_path = path.join(checkpointfolder, ui.replace("u", "w") + '.xdmf')
+            with XDMFFile(MPI.comm_world, checkpoint_path) as f:
+                f.write_checkpoint(w_[ui], '/current')
+
+            checkpoint_path = path.join(checkpointfolder, ui.replace("u", "d") + '.xdmf')
+            with XDMFFile(MPI.comm_world, checkpoint_path) as f:
+                f.write_checkpoint(d_[ui], '/current')
+            MPI.barrier(MPI.comm_world)
+
     # Store mesh and boundary
     MPI.barrier(MPI.comm_world)
     mesh_path = path.join(checkpointfolder, 'mesh.h5')
@@ -202,7 +216,7 @@ def check_if_reset_statistics(folder):
         return False
 
 
-def init_from_restart(restart_folder, sys_comp, uc_comp, u_components, q_, q_1, q_2, tstep, velocity_degree,
+def init_from_restart(restart_folder, sys_comp, uc_comp, u_components, q_, q_1, q_2, w_, d_, tstep, velocity_degree,
                       previous_velocity_degree, mesh, constrained_domain, V, Q, **NS_namespace):
     """Initialize solution from checkpoint files """
     if restart_folder:
@@ -219,6 +233,7 @@ def init_from_restart(restart_folder, sys_comp, uc_comp, u_components, q_, q_1,
             q_prev = dict((ui, Function(VV_prev[ui], name=ui)) for ui in sys_comp)
             q_2_prev = dict((ui, Function(V_prev, name=ui + "_2")) for ui in u_components)
 
+        # Load velocity and pressure
         for ui in sys_comp:
             checkpoint_path = path.join(restart_folder, ui + '.xdmf')
             with XDMFFile(MPI.comm_world, checkpoint_path) as f:
@@ -233,6 +248,16 @@ def init_from_restart(restart_folder, sys_comp, uc_comp, u_components, q_, q_1,
                         # Interpolate
                         read_and_interpolate_solution(f, V, previous_velocity_degree, q_2, q_2_prev, ui,
                                                       velocity_degree, "/previous")
+        # Load mesh velocity and deformation
+        if w_ is not None:
+            for ui in u_components:
+                checkpoint_path = path.join(restart_folder, ui.replace("u", "w") + '.xdmf')
+                with XDMFFile(MPI.comm_world, checkpoint_path) as f:
+                    f.read_checkpoint(w_[ui], "/current")
+
+                checkpoint_path = path.join(restart_folder, ui.replace("u", "d") + '.xdmf')
+                with XDMFFile(MPI.comm_world, checkpoint_path) as f:
+                    f.read_checkpoint(d_[ui], "/current")
 
 
 def read_and_interpolate_solution(f, V, previous_velocity_degree, q_, q_prev, ui, velocity_degree, name):

src/oasismove/solvers/NSfracStep/IPCS_ABCN_Move.py

@@ -283,7 +283,6 @@ def mesh_velocity_assemble(A_mesh, ui, bw, bw_tmp, a_mesh, bc_mesh, A_cache, **N
     A_mesh.zero()
     A_mesh.axpy(1, A_cache[(a_mesh, tuple(bc_mesh[ui]))], True)
     bw[ui].zero()
-    bw[ui].axpy(1., bw_tmp[ui])
 
 
 def mesh_velocity_solve(A_mesh, bw, wx_, w_, dof_map, dt, dx_, coordinates, w_sol, ui, bc_mesh, **NS_namespace):