Add lift/drag hooks, solution storage, plotting, and tests

pySDC/projects/StroemungsRaum/hooks/hooks_NSE_IMEX_FEniCS.py

@@ -0,0 +1,127 @@
+from pySDC.core.hooks import Hooks
+import dolfin as df
+import os
+
+
+class LogLiftDrag(Hooks):
+    """
+    Hook class to compute lift and drag forces on the cylinder at the end of each step,
+    and add them to the statistics
+    """
+
+    def post_step(self, step, level_number):
+        """
+        Record lift and drag forces on the cylinder at the end of each step
+
+        Args:
+            step (pySDC.Step.step): the current step
+            level_number (int): the current level number
+
+        Returns:
+            None
+        """
+        super().post_step(step, level_number)
+
+        L = step.levels[level_number]
+        L.sweep.compute_end_point()
+        P = L.prob
+
+        rho = 1
+
+        # normal pointing out of obstacle
+        n = -df.FacetNormal(P.V.mesh())
+
+        # tangential velocity component at the interface of the obstacle
+        u_t = df.inner(df.as_vector((n[1], -n[0])), L.uend.values)
+
+        # compute the drag and lift coefficients
+        drag = df.Form(2 / 0.1 * (P.nu / rho * df.inner(df.grad(u_t), n) * n[1] - P.pn * n[0]) * P.dsc)
+        lift = df.Form(-2 / 0.1 * (P.nu / rho * df.inner(df.grad(u_t), n) * n[0] + P.pn * n[1]) * P.dsc)
+
+        # assemble the scalar values
+        CD = df.assemble(drag)
+        CL = df.assemble(lift)
+
+        # add to statistics
+        self.add_to_stats(
+            process=step.status.slot,
+            time=L.time,
+            level=L.level_index,
+            iter=step.status.iter,
+            sweep=L.status.sweep,
+            type='lift_post_step',
+            value=CL,
+        )
+
+        self.add_to_stats(
+            process=step.status.slot,
+            time=L.time,
+            level=L.level_index,
+            iter=step.status.iter,
+            sweep=L.status.sweep,
+            type='drag_post_step',
+            value=CD,
+        )
+
+
+class LogWriteSolutions(Hooks):
+    """
+    Hook class to log the velocity and prssure solutions and write them in XDMF files
+    """
+
+    def pre_run(self, step, level_number):
+        """
+        Default routine called before time-loop starts
+
+        Args:
+            step (pySDC.Step.step): the current step
+            level_number (int): the current level number
+        """
+        super().pre_run(step, level_number)
+
+        L = step.levels[level_number]
+        P = L.prob
+
+        # create XDMF file for visualization output
+        path = f"{os.path.dirname(__file__)}/../data/navier_stokes/"
+        P.xdmffile_p = df.XDMFFile(path + 'Cylinder_pressure.xdmf')
+        P.xdmffile_u = df.XDMFFile(path + 'Cylinder_velocity.xdmf')
+
+    def post_step(self, step, level_number):
+        """
+        Write the solutions after each time step
+
+        Args:
+            step (pySDC.Step.step): the current step
+            level_number (int): the current level number
+
+        Returns:
+            None
+        """
+
+        super().post_step(step, level_number)
+
+        L = step.levels[level_number]
+        L.sweep.compute_end_point()
+        P = L.prob
+        t = L.time
+
+        P.xdmffile_p.write_checkpoint(P.pn, "pn", t, df.XDMFFile.Encoding.HDF5, True)
+        P.xdmffile_u.write_checkpoint(L.uend.values, "un", t, df.XDMFFile.Encoding.HDF5, True)
+
+    def post_run(self, step, level_number):
+        """
+        Default routine called after each run
+
+        Args:
+            step (pySDC.Step.step): the current step
+            level_number (int): the current level number
+        """
+
+        super().post_run(step, level_number)
+
+        L = step.levels[level_number]
+        P = L.prob
+
+        P.xdmffile_p.close()
+        P.xdmffile_u.close()

pySDC/projects/StroemungsRaum/plotting/plots_Navier_Stokes_projection_FEniCS.py

@@ -0,0 +1,135 @@
+import os
+import dolfin as df
+import json
+import matplotlib.animation as animation
+import matplotlib.pyplot as plt
+
+
+def plot_solutions():  # pragma: no cover
+
+    # Open XDMF files for reading the stored solutions
+    path = f"{os.path.dirname(__file__)}/../data/navier_stokes/"
+    xdmffile_u = df.XDMFFile(path + 'Cylinder_velocity.xdmf')
+    xdmffile_p = df.XDMFFile(path + 'Cylinder_pressure.xdmf')
+
+    # load parameters
+    parameters = json.load(open(path + "Navier_Stokes_FEniCS_parameters.json", 'r'))
+
+    # get parameters
+    dt = parameters['dt']
+    Tend = parameters['Tend']
+    family = parameters['family']
+    order = parameters['order']
+    t0 = parameters['t0']
+    nsteps = int(Tend / dt)
+
+    # load mesh
+    meshpath = f"{os.path.dirname(__file__)}/../meshs/"
+    mesh = df.Mesh(meshpath + 'cylinder.xml')
+
+    # define function spaces for velocity and physical quantities
+    V = df.VectorFunctionSpace(mesh, family, order)
+    Q = df.FunctionSpace(mesh, family, order - 1)
+    Vs = df.FunctionSpace(mesh, family, order)
+
+    # define functions for velocity and pressure
+    un = df.Function(V)
+    pn = df.Function(Q)
+
+    # stiffness  matrix for streamlines computation
+    u = df.TrialFunction(Vs)
+    v = df.TestFunction(Vs)
+    a = df.dot(df.nabla_grad(u), df.nabla_grad(v)) * df.dx
+    S = df.assemble(a)
+    Str = df.Function(Vs)
+
+    # open figure for plots
+    fig = plt.figure(1, figsize=(16, 13))
+
+    # initialize time variable
+    t = t0
+
+    for s in range(nsteps):
+
+        # Update current time
+        t += dt
+
+        # read the velocity and pressure solutions at the current time step
+        xdmffile_u.read_checkpoint(un, 'un', s)
+        xdmffile_p.read_checkpoint(pn, 'pn', s)
+
+        # compute the vorticity field
+        ux, uy = un.split(deepcopy=True)
+        Vort = uy.dx(0) - ux.dx(1)
+
+        # compute the streamlines
+        l = Vort * v * df.dx
+        L = df.assemble(l)
+        df.solve(S, Str.vector(), L)
+
+        # compute the magnitude of the velocity field
+        u_magn = df.sqrt(df.dot(un, un))
+        u_magn = df.project(u_magn, Vs)
+
+        # subplot for velocity
+        ax = fig.add_subplot(511)
+        c = df.plot(un, cmap='jet')
+        plt.colorbar(c)
+        ax.set_xlabel('Distance x')
+        ax.set_ylabel('Distance y')
+        ax.set_title('Velocity field')
+        ax.set_xlim(-0.01, 2.22)
+        ax.set_ylim(-0.005, 0.415)
+        plt.draw()
+
+        # subplot for pressure
+        ax = fig.add_subplot(512)
+        c = df.plot(pn, mode='color', cmap='jet')
+        plt.colorbar(c)
+        ax.set_xlabel('Distance x')
+        ax.set_ylabel('Distance y')
+        ax.set_title('pressure field')
+        ax.set_xlim(0, 2.20)
+        ax.set_ylim(0, 0.41)
+        plt.draw()
+
+        # subplot for vorticity
+        ax = fig.add_subplot(513)
+        c = df.plot(Vort, mode='color', vmin=-30, vmax=30, cmap='jet')
+        plt.colorbar(c)
+        ax.set_xlabel('Distance x')
+        ax.set_ylabel('Distance y')
+        ax.set_title('Vorticity')
+        ax.set_xlim(0, 2.20)
+        ax.set_ylim(0, 0.41)
+        plt.draw()
+
+        # subplot for velocity magnitude
+        ax = fig.add_subplot(514)
+        c = df.plot(u_magn, mode='color', cmap='jet')
+        plt.colorbar(c)
+        ax.set_xlabel('Distance x')
+        ax.set_ylabel('Distance y')
+        ax.set_title('Magnitude')
+        ax.set_xlim(0, 2.20)
+        ax.set_ylim(0, 0.41)
+        plt.draw()
+
+        ax = fig.add_subplot(515)
+        c = df.plot(Str, mode='contour', levels=50)
+        plt.colorbar(c)
+        ax.set_xlabel('Distance x')
+        ax.set_ylabel('Distance y')
+        ax.set_title('Magnitude = Streamlines')
+        ax.set_xlim(0, 2.20)
+        ax.set_ylim(0, 0.41)
+        plt.draw()
+
+        plt.pause(0.01)
+        plt.clf()
+
+    plt.close(fig)
+
+
+if __name__ == '__main__':
+    plot_solutions()

pySDC/projects/StroemungsRaum/run_Navier_Stokes_equations_FEniCS.py

@@ -1,6 +1,11 @@
+from pathlib import Path
+import json
+
 from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
 from pySDC.projects.StroemungsRaum.problem_classes.NavierStokes_2D_FEniCS import fenics_NSE_2D_mass
 from pySDC.projects.StroemungsRaum.sweepers.imex_1st_order_mass_NSE import imex_1st_order_mass_NSE
+from pySDC.projects.StroemungsRaum.hooks.hooks_NSE_IMEX_FEniCS import LogLiftDrag
+from pySDC.helpers.stats_helper import get_sorted
 
 
 def setup(t0=0):
@@ -46,6 +51,7 @@ def setup(t0=0):
     # initialize controller parameters
     controller_params = dict()
     controller_params['logger_level'] = 20
+    controller_params['hook_class'] = [LogLiftDrag]
 
     # Fill description dictionary
     description = dict()
@@ -95,6 +101,51 @@ def run_simulation(description, controller_params, Tend):
     return P, stats, uend
 
 
+def run_postprocessing(description, stats):
+    """
+    Postprocess and store simulation results for visualization and analysis.
+
+    Args:
+        description: dict,
+            pySDC description containing problem parameters.
+        problem: Problem instance,
+            Problem instance containing the final solution and other problem-related information.
+        stats: dict,
+            collected runtime statistics,
+
+    Returns: None
+    """
+    # get the data directory
+    import os
+
+    # create directory for storing results
+    path = f"{os.path.dirname(__file__)}/data/navier_stokes/"
+    Path(path).mkdir(parents=True, exist_ok=True)
+
+    # extract lift and drag coefficients from the collected statistics
+    lift = get_sorted(stats, type='lift_post_step', sortby='time')
+    drag = get_sorted(stats, type='drag_post_step', sortby='time')
+
+    # extract timing information
+    timing = get_sorted(stats, type='timing_run', sortby='time')
+
+    # save parameters
+    parameters = description['problem_params']
+    parameters.update(description['level_params'])
+    parameters['Tend'] = lift[-1][0]
+    parameters['timing'] = timing[0][1]
+    with open(path + "Navier_Stokes_FEniCS_parameters.json", 'w') as f:
+        json.dump(parameters, f)
+
+    # save lift and drag coefficients
+    with open(path + "Lift_Drag_Coefficients.txt", 'w') as f:
+        for i in range(len(lift)):
+            out = '%1.16f  %1.16f  %1.16f' % (lift[i][0], drag[i][1], lift[i][1])
+            f.write(out + '\n')
+
+    return None
+
+
 if __name__ == "__main__":
 
     t0 = 3.125e-04
@@ -105,3 +156,6 @@ def run_simulation(description, controller_params, Tend):
 
     # run the simulation and get the problem, stats and final solution
     problem, stats, uend = run_simulation(description, controller_params, Tend)
+
+    # run postprocessing to save parameters and solutions for visualization
+    run_postprocessing(description, stats)