Create incompressible_viscous solver and lid-driven cavity test probl…

…em (#138)

Needed to change a few things in incompressible/simulation.py as well:

Neumann BC's on phi in the case of dirichlet BC's on v (and apply phi-specific BC's in the multigrid solves)
Allow user-defined BC's
Allow another method for the velocity update in the evolve loop
The latter two are then redefined in incompressible_viscous/simulation.py. The rest of the methods are inherited from the base class.

For the lid-driven cavity problem, the only free parameter is the viscosity. Since the velocity and length scales are =1, the Reynolds number is just 1/viscosity.

.gitignore

@@ -22,7 +22,10 @@ inputs.auto
 
 __pycache__/
 
+.DS_Store
+
 
 pyro/_version.py
 pyro.egg-info/
+pyro_hydro.egg-info/
 .eggs/

README.md

@@ -157,6 +157,9 @@ pyro provides the following solvers (all in 2-d):
     projection method for the incompressible equations of
     hydrodynamics.
 
+  - `incompressible_viscous`: an extension of the incompressible solver
+    including a diffusion term for viscosity.
+
   - `lm_atm`: a solver for the equations of low Mach number
     hydrodynamics for atmospheric flows.
 

docs/source/incompressible_viscous_basics.rst

@@ -0,0 +1,115 @@
+Incompressible viscous hydrodynamics solver
+===========================================
+
+pyro's incompressible viscous solver solves:
+
+.. math::
+
+   \frac{\partial U}{\partial t} + U \cdot \nabla U + \nabla p &= \nu \nabla^2 U \\
+   \nabla \cdot U &= 0
+
+This solver is based on pyro's incompressible solver, but modifies the
+velocity update step to take viscosity into account.
+
+The main parameters that affect this solver are:
+
+.. include:: incompressible_viscous_defaults.inc
+
+Examples
+--------
+
+shear
+^^^^^
+
+The same shear problem as in incompressible solver, here with viscosity
+added.
+
+.. prompt:: bash
+
+   pyro_sim.py incompressible_viscous shear inputs.shear
+
+.. image:: shear_viscous.png
+   :align: center
+
+Compare this with the inviscid result. Notice how the velocities have 
+diffused in all directions.
+
+cavity
+^^^^^^
+
+The lid-driven cavity is a well-known benchmark problem for hydro codes
+(see e.g. :cite:t:`ghia1982`, :cite:t:`Kuhlmann2019`). In a unit square box 
+with initially static fluid, motion is initiated by a "lid" at the top 
+boundary, moving to the right with unit velocity. The basic command is:
+
+.. prompt:: bash
+
+   pyro_sim.py incompressible_viscous cavity inputs.cavity
+
+It is interesting to observe what happens when varying the viscosity, or, 
+equivalently the Reynolds number (in this case :math:`\rm{Re}=1/\nu` since 
+the characteristic length and velocity scales are 1 by default).
+
+|pic1| |pic2| |pic3|
+
+.. |pic1| image:: cavity_Re100.png
+   :width: 32%
+
+.. |pic2| image:: cavity_Re400.png
+   :width: 32%
+
+.. |pic3| image:: cavity_Re1000.png
+   :width: 32%
+
+These plots were made by allowing the code to run for longer and approach a
+steady-state with the option ``driver.max_steps=1000``, then running 
+(e.g. for the Re=100 case):
+
+.. prompt:: bash
+
+   python incompressible_viscous/problems/plot_cavity.py cavity_n64_Re100_0406.h5 -Re 100 -o cavity_Re100.png
+
+convergence
+^^^^^^^^^^^
+
+This is the same test as in the incompressible solver. With viscosity,
+an exponential term is added to the solution. Limiting can again be 
+disabled by adding ``incompressible.limiter=0`` to the run command.
+The basic set of tests shown below are run as:
+
+.. prompt:: bash
+
+   pyro_sim.py incompressible_viscous converge inputs.converge.32 vis.dovis=0
+   pyro_sim.py incompressible_viscous converge inputs.converge.64 vis.dovis=0
+   pyro_sim.py incompressible_viscous converge inputs.converge.128 vis.dovis=0
+
+The error is measured by comparing with the analytic solution using
+the routine ``incomp_viscous_converge_error.py`` in ``analysis/``. To 
+generate the plot below, run
+
+.. prompt:: bash
+
+   python incompressible_viscous/tests/convergence_errors.py convergence_errors.txt
+
+or ``convergence_errors_no_limiter.txt`` after running with that option. Then:
+
+.. prompt:: bash
+
+   python incompressible_viscous/tests/convergence_plot.py
+
+.. image:: incomp_viscous_converge.png
+   :align: center
+
+The solver is converging but below second-order, unlike the inviscid case. Limiting
+does not seem to make a difference here.
+
+Exercises
+---------
+
+Explorations
+^^^^^^^^^^^^
+
+* In the lid-driven cavity problem, when does the solution reach a steady-state?
+
+* :cite:`ghia1982` give benchmark velocities at different Reynolds number for the
+  lid-driven cavity problem (see their Table I). Do we agree with their results?

docs/source/incompressible_viscous_defaults.inc

@@ -0,0 +1,36 @@
+* section: [driver]
+
+  +----------------------------------+----------------+----------------------------------------------------+
+  | option                           | value          | description                                        |
+  +==================================+================+====================================================+
+  | cfl                              | ``0.8``        |                                                    |
+  +----------------------------------+----------------+----------------------------------------------------+
+
+* section: [incompressible]
+
+  +----------------------------------+----------------+----------------------------------------------------+
+  | option                           | value          | description                                        |
+  +==================================+================+====================================================+
+  | limiter                          | ``2``          | limiter (0 = none, 1 = 2nd order, 2 = 4th order)   |
+  +----------------------------------+----------------+----------------------------------------------------+
+  | proj_type                        | ``2``          | what are we projecting? 1 includes -Gp term in U*  |
+  +----------------------------------+----------------+----------------------------------------------------+
+
+* section: [incompressible_viscous]
+
+  +----------------------------------+----------------+----------------------------------------------------+
+  | option                           | value          | description                                        |
+  +==================================+================+====================================================+
+  | viscosity                        | ``0.1``        | kinematic viscosity of the fluid (units L^2/T)     |
+  +----------------------------------+----------------+----------------------------------------------------+
+
+* section: [particles]
+
+  +----------------------------------+----------------+----------------------------------------------------+
+  | option                           | value          | description                                        |
+  +==================================+================+====================================================+
+  | do_particles                     | ``0``          |                                                    |
+  +----------------------------------+----------------+----------------------------------------------------+
+  | particle_generator               | ``grid``       |                                                    |
+  +----------------------------------+----------------+----------------------------------------------------+
+

docs/source/index.rst

@@ -36,6 +36,7 @@ pyro: a python hydro code
    multigrid
    diffusion_basics
    incompressible_basics
+   incompressible_viscous_basics
    lowmach_basics
    swe_basics
    particles_basics

docs/source/pyro.incompressible_viscous.problems.rst

@@ -0,0 +1,42 @@
+pyro.incompressible\_viscous.problems package
+=============================================
+
+.. automodule:: pyro.incompressible_viscous.problems
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+Submodules
+----------
+
+pyro.incompressible\_viscous.problems.cavity module
+---------------------------------------------------
+
+.. automodule:: pyro.incompressible_viscous.problems.cavity
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+pyro.incompressible\_viscous.problems.converge module
+-----------------------------------------------------
+
+.. automodule:: pyro.incompressible_viscous.problems.converge
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+pyro.incompressible\_viscous.problems.plot\_cavity module
+---------------------------------------------------------
+
+.. automodule:: pyro.incompressible_viscous.problems.plot_cavity
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+pyro.incompressible\_viscous.problems.shear module
+--------------------------------------------------
+
+.. automodule:: pyro.incompressible_viscous.problems.shear
+   :members:
+   :undoc-members:
+   :show-inheritance:

docs/source/pyro.incompressible_viscous.rst

@@ -0,0 +1,34 @@
+pyro.incompressible\_viscous package
+====================================
+
+.. automodule:: pyro.incompressible_viscous
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+Subpackages
+-----------
+
+.. toctree::
+   :maxdepth: 4
+
+   pyro.incompressible_viscous.problems
+
+Submodules
+----------
+
+pyro.incompressible\_viscous.BC module
+--------------------------------------
+
+.. automodule:: pyro.incompressible_viscous.BC
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+pyro.incompressible\_viscous.simulation module
+----------------------------------------------
+
+.. automodule:: pyro.incompressible_viscous.simulation
+   :members:
+   :undoc-members:
+   :show-inheritance:

docs/source/pyro.rst

@@ -26,6 +26,7 @@ Subpackages
    pyro.compressible_sr
    pyro.diffusion
    pyro.incompressible
+   pyro.incompressible_viscous
    pyro.lm_atm
    pyro.mesh
    pyro.multigrid

docs/source/refs.bib

@@ -37,7 +37,6 @@ @article{ZALESAK1979335
 abstract = "The theory of flux-corrected transport (FCT) developed by Boris and Book [J. Comput. Phys. 11 (1973) 38; 18 (1975) 248; 20 (1976) 397] is placed in a simple, generalized format, and a new algorithm for implementing the critical flux limiting stage' in multidimensions without resort to time splitting is presented. The new flux limiting algorithm allows the use of FCT techniques in multidimensional fluid problems for which time splitting would produce unacceptable numerical results, such as those involving incompressible or nearly incompressible flow fields. The “clipping” problem associated with the original one dimensional flux limiter is also eliminated or alleviated. Test results and applications to a two dimensional fluid plasma problem are presented."
 }
 
-
 @article{stoker:1958,
   author = {Stoker,J. J.  and Lindsay,R. Bruce },
   title = {Water Waves},
@@ -101,3 +100,32 @@ @ARTICLE{minion:1996
        adsurl = {https://ui.adsabs.harvard.edu/abs/1996JCoPh.127..158M},
       adsnote = {Provided by the SAO/NASA Astrophysics Data System}
 }
+
+@ARTICLE{ghia1982,
+       author = {{Ghia}, U. and {Ghia}, K.~N. and {Shin}, C.~T.},
+        title = "{High-Re Solutions for Incompressible Flow Using the Navier-Stokes Equations and a Multigrid Method}",
+      journal = {Journal of Computational Physics},
+     keywords = {Computational Fluid Dynamics, Computational Grids, High Reynolds Number, Incompressible Flow, Navier-Stokes Equation, Algorithms, Computer Programs, Partial Differential Equations, Velocity Distribution, Fluid Mechanics and Heat Transfer},
+         year = 1982,
+        month = dec,
+       volume = {48},
+       number = {3},
+        pages = {387-411},
+          doi = {10.1016/0021-9991(82)90058-4},
+       adsurl = {https://ui.adsabs.harvard.edu/abs/1982JCoPh..48..387G},
+      adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+@Inbook{Kuhlmann2019,
+  author="Kuhlmann, Hendrik C.
+  and Roman{\`o}, Francesco",
+  editor="Gelfgat, Alexander",
+  title="The Lid-Driven Cavity",
+  bookTitle="Computational Modelling of Bifurcations and Instabilities in Fluid Dynamics",
+  year="2019",
+  publisher="Springer International Publishing",
+  address="Cham",
+  pages="233--309",
+  isbn="978-3-319-91494-7",
+  doi="10.1007/978-3-319-91494-7_8",
+}

pyro/analysis/incomp_converge_error.py

@@ -47,12 +47,7 @@ def main():
         print(usage)
         sys.exit(2)
 
-    try:
-        file = sys.argv[1]
-    except IndexError:
-        print(usage)
-        sys.exit(2)
-
+    file = sys.argv[1]
     errors = get_errors(file)
     print("errors: ", errors[0], errors[1])
 

pyro/analysis/incomp_viscous_converge_error.py

@@ -0,0 +1,57 @@
+#!/usr/bin/env python3
+
+
+import math
+import sys
+
+import numpy as np
+
+import pyro.util.io_pyro as io
+
+usage = """
+      compare the output in file from the incompressible viscous converge problem to
+      the analytic solution.
+
+      usage: ./incomp_viscous_converge_error.py file
+"""
+
+
+def get_errors(file):
+
+    sim = io.read(file)
+    myd = sim.cc_data
+    myg = myd.grid
+
+    # numerical solution
+    u = myd.get_var("x-velocity")
+    v = myd.get_var("y-velocity")
+    nu = myd.get_aux("viscosity")
+
+    t = myd.t
+
+    # analytic solution
+    u_exact = myg.scratch_array()
+    u_exact[:, :] = 1.0 - 2.0*np.cos(2.0*math.pi*(myg.x2d-t))*np.sin(2.0*math.pi*(myg.y2d-t))*np.exp(-8*math.pi**2*nu*t)
+
+    v_exact = myg.scratch_array()
+    v_exact[:, :] = 1.0 + 2.0*np.sin(2.0*math.pi*(myg.x2d-t))*np.cos(2.0*math.pi*(myg.y2d-t))*np.exp(-8*math.pi**2*nu*t)
+
+    # error
+    udiff = u_exact - u
+    vdiff = v_exact - v
+
+    return udiff.norm(), vdiff.norm()
+
+
+def main():
+    if len(sys.argv) != 2:
+        print(usage)
+        sys.exit(2)
+
+    file = sys.argv[1]
+    errors = get_errors(file)
+    print("errors: ", errors[0], errors[1])
+
+
+if __name__ == "__main__":
+    main()