Code cleaning and math fix (KVSlab#2)

* Fix missing + and move theta outside of Piola. * Code style clean, and add comments. * Add comments to the code. * Fix typo * More fancy verbose=False print * Update problem files * Update print function and move it to __init__ for further user flexebility * Update reference values * Fix formating error in docs * Update restart_folder check * Update documentation. Note: Missing HDF5 parallel fix. * visualization and xdmf fix * doc update
alexruiz95 · Dec 5, 2019 · 04faa47 · 04faa47
1 parent c053d31
commit 04faa47
Show file tree

Hide file tree

Showing 19 changed files with 363 additions and 345 deletions.
diff --git a/docs/source/acknow_ref.rst b/docs/source/acknow_ref.rst
@@ -13,6 +13,6 @@ The numerical schemes in turtleFSI is presented and tested in Slyngstad [1]_ and
 The input problem set up for the TF_cfd, TF_csm, and TF_fsi is taken from the Turek et al. [3]_ benchmark
 paper.
 
-.. [1] Slyngstad, Andreas Strøm. Verification and Validation of a Monolithic Fluid-Structure Interaction Solver in FEniCS. A comparison of mesh lifting operators. MS thesis. 2017.
+.. [1] Slyngstad, Andreas S. Verification and Validation of a Monolithic Fluid-Structure Interaction Solver in FEniCS. A comparison of mesh lifting operators. MS thesis. 2017.
 .. [2] Gjertsen, Sebastian. Development of a Verified and Validated Computational Framework for Fluid-Structure Interaction: Investigating Lifting Operators and Numerical Stability. MS thesis. 2017.
 .. [3] Turek, Stefan, and Jaroslav Hron. "Proposal for numerical benchmarking of fluid-structure interaction between an elastic object and laminar incompressible flow." Fluid-structure interaction. Springer, Berlin, Heidelberg, 2006. 371-385.
diff --git a/docs/source/conf.py b/docs/source/conf.py
@@ -24,9 +24,9 @@
 author = 'A. Sylngstad, S. Gjertsen, A. Bergersen, A. Souche, and K. Valen-Sendstad'
 
 # The short X.Y version
-version = 'v1.1'
+version = 'v1.2'
 # The full version, including alpha/beta/rc tags
-release = 'v1.1.0'
+release = 'v1.2.0'
 
 
 # -- General configuration ---------------------------------------------------

diff --git a/docs/source/installation.rst b/docs/source/installation.rst
@@ -8,7 +8,7 @@ Installation
 
 Compatibility and Dependencies
 ==============================
-The general dependencies of turtleFSI are:
+The dependencies of turtleFSI are:
 
 * FEniCS 2019.1.0
 * Numpy >1.1X
@@ -23,14 +23,15 @@ Then execute the following command in a terminal window::
     $ conda create -n your_environment -c conda-forge turtleFSI
 
 You can then activate your environment by running ``source activate your_environment``.
-Now you are all set, and can start using turtleFSI. You can execute turtleFSI directly
-through the terminal by typing::
+Now you are all set, and can start using turtleFSI. A detailed explanation for usage of
+turtleFSI can be found here_.
 
-    $ turtleFSI
+If you are using turtleFSI on a high performance computing (HPC) cluster we always
+recommend that you build from source, as described below. This is in accordance
+with the guidelines provided by the `FEniCS project <https://fenicsproject.org/download/>`_
+users to install FEniCS from source when on a HPC cluster.
 
-followed by any additional options, for instance, which problem to run and the time step size.
-Use ``-h`` to see all available options. A detailed explanation for usage of turtleFSI can be found
-`here <https://turtlefsi2.readthedocs.io/en/latest/using_turtleFSI.html>`.
+.. _here: <https://turtlefsi2.readthedocs.io/en/latest/using_turtleFSI.html>.
 
 
 Development version
@@ -40,12 +41,12 @@ Downloading
 ~~~~~~~~~~~
 The latest development version of turtleFSI can be found on the official
 `turtleFSI git repository <https://github.com/KVSlab/turtleFSI>`_ on Github.
-To clone the turtleFSI repository, navigate to the directory where you wish
+To clone the turtleFSI repository, open a terminal, navigate to the directory where you wish
 turtleFSI to be stored, type the following command, and press Enter::
 
     $ git clone https://github.com/KVSlab/turtleFSI
 
-After the source distribution has been downloaded, all the files required will be located
+After the source distribution has been downloaded, all the files will be located
 in the newly created ``turtleFSI`` folder.
 
 Building
@@ -55,3 +56,7 @@ file will be located. First, make sure that all dependencies are installed.
 Then, you can install turtleFSI be executing the following::
 
     $ python setup.py install
+
+
+If you are installing turtleFSI somewhere you do not have root access, typically on a cluster, you can add
+``--user`` to install locally.
diff --git a/docs/source/using_turtleFSI.rst b/docs/source/using_turtleFSI.rst
@@ -11,37 +11,47 @@ Execute a run
 =============
 
 turtleFSI is aimed to be user friendly, where all parameters can be controlled from the command line.
-In this page, we provide an example on how to create your own problem file. First a quick recap on how to
+We here provide an example on how to create your own problem file. First a quick recap on how to
 execute turtleFSI.
 
 To run turtleFSI with all the default parameters you can execute::
+
  turtleFSI
 
-or, to run a specific problem file, run::
+which will execute the default problem file ``turtle_demo.py`` with all the default parameter. In the folder where
+you executed the command there will be a folder ``turtle_demo_results/1`` where you can find ``Visualization`` files,
+and under ``Checkpoint`` you can find a list of all parameters used in ``default_variables.pickle``.
+
+To run a specific problem file, run::
+
  turtleFSI --problem [path_to_problem]
 
 To get an overview of all parameters, please run::
+
  turtleFSI -h
 
 
 Built-in functionality
 ======================
-Two important functions are available in the current version of turtleFSI: checkpointing/restart, and storing files for visualization.
+TurtleFSI is designed to be lightweight and generic, but we have added to features generally of interest: checkpointing/restart, and
+storing files for visualization.
 
-For checkpointing you can set the variable ``--checkpoint-step`` to set how often a checkpoint should
-be stored. To restart from a previous checkpoint step use the command ``--restart-folder [folder/sub-folder]``. Note that
-the variables from the previous simulation will overwrite any parameters set in the ``set_problem_parameters``
-or on the commandline. If you need to change a parameter from the previous checkpoint file (for instance, end time ``T``), you can
-still do it by explicitly redefining it within the ``initiate`` function.
+For checkpointing you can set the variable ``--checkpoint-step`` to set the checkpoint frequency, i.e., how often a
+checkpoint should be stored. To restart from a previous checkpoint, use the command
+``--restart-folder [folder]/[sub-folder]``. Note that the variables from the previous simulation will overwrite any
+parameters defined in ``set_problem_parameters`` or on the commandline. If you need to change a parameter from the
+previous checkpoint file (for instance, end time ``T``), you can still do it by explicitly redefining the variable
+in the ``initiate`` function.
 
-To set how often you save files for visualization you can set ``--save-step``. Note that the default is ``1``.
+To set how often you save files for visualization you can set ``--save-step``. Note that the default is ``10``.
 
 
 Setting parameters
 ==================
 All the default parameters are set in the ``problem/__init__.py`` file. Problem specific parameters
-are then overwritten in the problem file under ``set_problem_parameters`` or by defining them in the command line. In summary;
-the priority is as follow: default parameters < problem file < command line < (checkpointing).
+are then overwritten in the problem file under ``set_problem_parameters`` or from the command line.
+In summary; the parameters will be updated in the following order: default parameters < ``set_problem_parameters`` <
+command line (< checkpointing) < other problem specific functions.
 
 
 Create your own problem file
@@ -56,27 +66,23 @@ the above mentioned information. Listed in the order they are first executed:
 
 - ``set_problem_parameters``
 - ``get_mesh_domain_and_boundaries``
-- ``initiate``
+- ``initiate`` (optional)
 - ``create_bcs``
-- ``pre_solve``
-- ``post_solve``
-- ``finished``
+- ``pre_solve`` (optional)
+- ``post_solve`` (optional)
+- ``finished`` (optional)
 
 
 set_problem_parameters
 ~~~~~~~~~~~~~~~~~~~~~~
-This function is for defining parameters of the problem like, dt, end time (T), and
-physical parameters of the problem. To see a full list of the standard parameters you can change
-please refer to the ``default_variables`` defined in ``turtleFSI/problems/__init__.py``.
+This function is for defining parameters of the problem like: time step size ``dt``, end time ``T``, and
+physical parameters of the problem. To see a full list of the default parameters, we refer to the
+``default_variables`` defined in ``turtleFSI/problems/__init__.py``. Please note that you can, and should, also
+define other problem specific variables in ``default_variables``, for instance, geometry information like hight
+and length of the problem is of interest to the problem you are solving.
 
-In ``set_problem_parameters`` you should take the ``default_variables`` as an input,
-and update the dictionary with your own values. It is particularly important to
-overwrite the physical variables as these vary from problem to problem.
-
-If you provide any command line arguments these will overwrite both those you have defined in your
-problem file, and the ``default_variables``. In theory, you do not have to specify the ``set_problem_parameters``
-if you just want to use the values defined in ``default_variables``, however in practice you have to
-include this function in all your problem files.
+In the ``set_problem_parameters`` function, you have the possibility to set your problem parameters within the
+``default_variables`` dictionary. However, keep in mind that any command line arguments will overwrite the ``default_variables``.
 
 A simple example of this function can look like this::
 
@@ -104,15 +110,21 @@ A simple example of this function can look like this::
          return default_variables
 
 
+.. note::
+    The parameter extrapolation here refers to mesh lifting operator, i.e., how to de extrapolate the deformation of
+    the solid into the fluid to create an ALE frame of reference. Laplace is the cheapest in terms of computational
+    cost, but is less robust for large deformations and sharp edges. In contrast, biharmonic is very robust, but
+    at the cost of computational efficiency and increased memory load.
+
 get_mesh_domain_and_boundaries
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-This function is one of two which is not optional, and have to be provided for the problem file to run.
-Here you read or define your mesh, domain markers, and boundary markers. In ``turtle_demo.py`` there
-is an example of reading the mesh data from pre-existing ".xdmf" mesh files. In contrast to the
-'Turek flag'-example (``TF_fsi.py``), where the domain and boundaries are marked using FEniCS functions.
+The function is essential and unique for every problem, and has to be provided for the problem file to run.
+Here you read or define your mesh, domain markers, and boundary markers. In ``turtle_demo.py`` we
+read the mesh data from pre-existing ".xdmf" mesh files. In contrast to defining the domain and
+boundary markers using FEniCS functions, like in the the 'Turek flag'-example (``TF_fsi.py``).
 
-Any questions regarding how to best create a mesh, we refer to the FEniCS documentation and discourse group, but
-``pygmsh`` in combination with ``meshio`` can be relevant tools to create a lot of geometries.
+``pygmsh`` and ``meshio`` are relevant tools to create geometries. For any question regarding meshing,
+we refer to the FEniCS documentation and `discourse group <https://fenicsproject.discourse.group/>`_.
 
 
 In ``turtle_demo.py``, the function looks like this::
@@ -172,7 +184,6 @@ This class is then used in the function ``create_bcs`` to prescribe Dirichlet bo
 inlet velocity. When defining the boundary conditions to specific domain regions or boundaries, make sure to
 be consistent with the markers provided in ``get_mesh_domain_and_boundaries``::
 
-
     class Inlet(UserExpression):
         def __init__(self, Um, **kwargs):
             self.t = 0.0
@@ -214,8 +225,8 @@ be consistent with the markers provided in ``get_mesh_domain_and_boundaries``::
 
         # Fluid velocity boundary conditions
         u_inlet = DirichletBC(DVP.sub(1), inlet, boundaries, inlet_id)
-        u_bot = DirichletBC(DVP.sub(1).sub(1), (0.0), boundaries, bottom_id)
-        u_top = DirichletBC(DVP.sub(1).sub(1), (0.0), boundaries, top_id)
+        u_bot = DirichletBC(DVP.sub(1).sub(1), (0.0), boundaries, bottom_id)  # slip in x-direction
+        u_top = DirichletBC(DVP.sub(1).sub(1), (0.0), boundaries, top_id)     # slip in x-direction
         u_head_tail = DirichletBC(DVP.sub(1), noslip, boundaries, turtle_head_tail_id)
 
         # Pressure boundary conditions
@@ -236,12 +247,11 @@ be consistent with the markers provided in ``get_mesh_domain_and_boundaries``::
 
         return dict(bcs=bcs, inlet=inlet)
 
-
 .. figure:: ../../figs/Turtle_boundaries_zoom.png
     :width: 600px
     :align: center
 
-    FSI and Fixed boundaries.
+    Boundaries between the fluid and structures and fixed boundaries.
 
 .. figure:: ../../figs/Turtle_inlet_vel.png
    :width: 600px

diff --git a/tests/test_turtleFSI.py b/tests/test_turtleFSI.py
@@ -23,8 +23,8 @@ def test_cfd():
 
     drag = np.loadtxt(Path.cwd().joinpath("tmp/1/Drag.txt"))[-1]
     lift = np.loadtxt(Path.cwd().joinpath("tmp/1/Lift.txt"))[-1]
-    drag_reference = 2.5637554331614054
-    lift_reference = -0.02078995609237899
+    drag_reference = 4.503203576965564
+    lift_reference = -0.03790359084395478
 
     assert compare(drag, drag_reference)
     assert compare(lift, lift_reference)
@@ -37,8 +37,8 @@ def test_csm():
 
     distance_x = np.loadtxt("tmp/2/dis_x.txt")[-1]
     distance_y = np.loadtxt("tmp/2/dis_y.txt")[-1]
-    distance_x_reference = -6.13487990897633e-06
-    distance_y_reference = -3.9398599897576816e-05
+    distance_x_reference = -3.312418050495862e-05
+    distance_y_reference = -0.003738529237136441
 
     assert compare(distance_x, distance_x_reference)
     assert compare(distance_y, distance_y_reference)
@@ -53,10 +53,10 @@ def test_fsi():
     lift = np.loadtxt("tmp/3/Lift.txt")[-1]
     distance_x = np.loadtxt("tmp/3/dis_x.txt")[-1]
     distance_y = np.loadtxt("tmp/3/dis_y.txt")[-1]
-    distance_x_reference = -3.0193475393178104e-06
-    distance_y_reference = -2.6594621765280973e-08
-    drag_reference = 2.472928697982334
-    lift_reference = -0.003950732063025264
+    distance_x_reference = -6.896013956339182e-06
+    distance_y_reference = 1.876355330341896e-09
+    drag_reference = 4.407481239804155
+    lift_reference = -0.005404703556977697
 
     assert compare(distance_x, distance_x_reference)
     assert compare(distance_y, distance_y_reference)
@@ -75,10 +75,10 @@ def test_laplace(extrapolation_sub_type):
     lift = np.loadtxt("tmp/4/Lift.txt")[-1]
     distance_x = np.loadtxt("tmp/4/dis_x.txt")[-1]
     distance_y = np.loadtxt("tmp/4/dis_y.txt")[-1]
-    distance_x_reference = -3.0193475393178104e-06
-    distance_y_reference = -2.6594621765280973e-08
-    drag_reference = 2.472928697982334
-    lift_reference = -0.003950732063025264
+    distance_x_reference = -6.896013956339182e-06
+    distance_y_reference = 1.876355330341896e-09
+    drag_reference = 4.407481239804155
+    lift_reference = -0.005404703556977697
 
     assert compare(distance_x, distance_x_reference)
     assert compare(distance_y, distance_y_reference)
@@ -96,10 +96,10 @@ def test_biharmonic(extrapolation_sub_type):
     lift = np.loadtxt("tmp/5/Lift.txt")[-1]
     distance_x = np.loadtxt("tmp/5/dis_x.txt")[-1]
     distance_y = np.loadtxt("tmp/5/dis_y.txt")[-1]
-    distance_x_reference = -3.0193475393178104e-06
-    distance_y_reference = -2.6594621765280973e-08
-    drag_reference = 2.472928697982334
-    lift_reference = -0.003950732063025264
+    distance_x_reference = -6.896013956339182e-06
+    distance_y_reference = 1.876355330341896e-09
+    drag_reference = 4.407481239804155
+    lift_reference = -0.005404703556977697
 
     assert compare(distance_x, distance_x_reference)
     assert compare(distance_y, distance_y_reference)
@@ -116,10 +116,10 @@ def test_elastic():
     lift = np.loadtxt("tmp/6/Lift.txt")[-1]
     distance_x = np.loadtxt("tmp/6/dis_x.txt")[-1]
     distance_y = np.loadtxt("tmp/6/dis_y.txt")[-1]
-    distance_x_reference = -3.019356900018008e-06
-    distance_y_reference = -2.659921700576888e-08
-    drag_reference = 2.472926160030297
-    lift_reference = -0.0039508769183140835
+    distance_x_reference = -6.896144755254494e-06
+    distance_y_reference = 1.868651990487361e-09
+    drag_reference = 4.407488867909029
+    lift_reference = -0.005404616050528832
 
     assert compare(distance_x, distance_x_reference)
     assert compare(distance_y, distance_y_reference)

diff --git a/turtleFSI/modules/biharmonic.py b/turtleFSI/modules/biharmonic.py
@@ -3,19 +3,19 @@
 # the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
 # PURPOSE.
 
-from dolfin import inner, grad, div
+from dolfin import inner, grad
 
 
 def extrapolate_setup(F_fluid_linear, extrapolation_sub_type, d_, w_, phi, beta, dx_f,
                       ds, n, bc_ids, **namespace):
     """
     Biharmonic lifting operator. Should be used for large deformations.
 
-    alfa * laplace^2(d) = 0   in the fluid domain
+    alpha * laplace^2(d) = 0   in the fluid domain
 
     By introducing w = - grad(d) we obtain the equivalent system of equations:
-        w = - alfa * laplace(d)
-        - alfa * grad(w) = 0
+        w = - alpha * laplace(d)
+        - alpha * grad(w) = 0
 
     Two types of boundary conditions can be setup for this problem:
 
@@ -26,27 +26,16 @@ def extrapolate_setup(F_fluid_linear, extrapolation_sub_type, d_, w_, phi, beta,
         - "constrained_disp_vel" with conditions on (d) and (w):
             d(d(x))/dn = 0 and d(w(x))/dn = 0   on the inlet and outlet fluid boundaries
             d(d(y))/dn = 0 and d(w(y))/dn = 0   on the FSI interface
-
-    References:
-
-    Slyngstad, Andreas Strøm. Verification and Validation of a Monolithic
-        Fluid-Structure Interaction Solver in FEniCS. A comparison of mesh lifting
-        operators. MS thesis. 2017.
-
-    Gjertsen, Sebastian. Development of a Verified and Validated Computational
-        Framework for Fluid-Structure Interaction: Investigating Lifting Operators
-        and Numerical Stability. MS thesis. 2017.
     """
 
-    alfa_u = 0.01
-    F_ext1 = alfa_u*inner(w_["n"], beta)*dx_f - alfa_u*inner(grad(d_["n"]),
-                                                             grad(beta))*dx_f
-    F_ext2 = alfa_u*inner(grad(w_["n"]), grad(phi))*dx_f
+    alpha_u = 0.01
+    F_ext1 = alpha_u * inner(w_["n"], beta) * dx_f - alpha_u * inner(grad(d_["n"]), grad(beta)) * dx_f
+    F_ext2 = alpha_u * inner(grad(w_["n"]), grad(phi)) * dx_f
 
     if extrapolation_sub_type == "constrained_disp_vel":
         for i in bc_ids:
-            F_ext1 += alfa_u*inner(grad(d_["n"])*n, beta)*ds(i)
-            F_ext2 -= alfa_u*inner(grad(w_["n"])*n, phi)*ds(i)
+            F_ext1 += alpha_u*inner(grad(d_["n"]) * n, beta) * ds(i)
+            F_ext2 -= alpha_u*inner(grad(w_["n"]) * n, phi) * ds(i)
 
     F_fluid_linear += F_ext1 + F_ext2
 

diff --git a/turtleFSI/modules/elastic.py b/turtleFSI/modules/elastic.py
@@ -14,22 +14,12 @@ def extrapolate_setup(F_fluid_linear, mesh, d_, phi, gamma, dx_f, **namespace):
     div(sigma(d)) = 0   in the fluid domain
     d = 0               on the fluid boundaries other than FSI interface
     d = solid_d         on the FSI interface
-
-    References:
-
-    Slyngstad, Andreas Strøm. Verification and Validation of a Monolithic
-        Fluid-Structure Interaction Solver in FEniCS. A comparison of mesh lifting
-        operators. MS thesis. 2017.
-
-    Gjertsen, Sebastian. Development of a Verified and Validated Computational
-        Framework for Fluid-Structure Interaction: Investigating Lifting Operators
-        and Numerical Stability. MS thesis. 2017.
     """
-    E_y = 1./CellVolume(mesh)
+    E_y = 1.0 / CellVolume(mesh)
     nu = 0.25
-    alfa_lam = nu*E_y / ((1. + nu)*(1. - 2.*nu))
-    alfa_mu = E_y/(2.*(1. + nu))
-    F_extrapolate = inner(J_(d_["n"]) * S_linear(d_["n"], alfa_mu, alfa_lam) *
+    alpha_lam = nu * E_y / ((1.0 + nu) * (1.0 - 2.0 * nu))
+    alpha_mu = E_y / (2.0 * (1.0 + nu))
+    F_extrapolate = inner(J_(d_["n"]) * S_linear(d_["n"], alpha_mu, alpha_lam) *
                           inv(F_(d_["n"])).T, grad(phi))*dx_f
 
     F_fluid_linear += F_extrapolate