First commit, included main, solvers, and ipynb

.gitignore

@@ -0,0 +1,41 @@
+# Byte-compiled / optimized / DLL files
+__pycache__/
+
+# PEP 582; used by e.g. github.com/David-OConnor/pyflow
+__pypackages__/
+
+# Distribution / packaging
+.Python
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+wheels/
+share/python-wheels/
+*.egg-info/
+.installed.cfg
+*.egg
+MANIFEST
+
+# Generated results from FEniCS/Python plots, etc...
+*.xdmf
+*.h5
+*.pyc
+*.png
+
+# Jupyter checkpoints
+#*.ipynb
+.ipynb_checkpoints/
+
+# Numerical results
+tests/
+
+# Paraview state files
+*.pvsm

README.rst

@@ -0,0 +1,63 @@
+NavierStokes
+============
+
+This repository implements via finite element solvers for incompressible
+Navier-Stokes (iNSE) equations in Arbitrary Lagrangian-Eulerian (ALE) formalism the schemes 
+proposed on [RA20]_.
+
+In particular, it aims to replicate the energy results shown in the article [RA20]_ 
+for both Monolithic and Chorin-Temam solvers.
+
+Solvers
+-------
+
+Both monolithic solver and a fractional step solver are implemented.
+
+* Monolithic solver for the iNSE-ALE problem with linearized convective term
+and Taylor-Hood (P2/P1) stable finite element space.
+
+
+* Fractional step solver for the iNSE-ALE problem with linealized convective term
+and P1/P1 finite element space. The Chorin-Temam schemes proposed is described in [RA20]_.
+
+Flow model
+----------
+
+A rectangle domain is taken with fully Dirichlet homogeneous boundary conditions
+and non-zero initial velocity profile. Further description of the problem also can be found 
+in the reference.
+
+Usage
+----------
+
+Since the repository aims to directly reproduce the results of the reference, 
+no configuration files where implemented to further customize the problem. 
+Nevertheless, the solvers are easily modified since its implementation is done via 
+FEniCS [LO12]_.
+
+An ipynb file is included to reproduce the results and recommended to use.
+
+If desired to run the simulations only, execute::
+
+    python main.py
+
+
+Dependencies
+------------
+
+- Python >= 3.5 
+- FEniCS_ >= 2019.1.0
+
+
+
+Reference
+^^^^^^^^^^
+
+.. [RA20] Aróstica R., Bertoglio C. (2020) On monolithic and Chorin-Temam
+   schemes for incompressible flows in moving domains. 
+   Applied Mathematics Letters, doi: https://doi.org/10.1016/j.aml.2020.106830
+   ISBN: 978-3-642-23099-8
+
+.. [LO12] Logg A., Mardal K.-A., Wells G. N. (2012) Automated Solution of Differential
+   Equations by the Finite Element Method.
+   Springer, Berlin, Heidelberg, doi: https://doi.org/10.1007/978-3-642-23099-8

common.py

@@ -0,0 +1,54 @@
+"""
+This module contains all common functions used in variational form. 
+All operators are defined with respect to the reference coordiantes.
+"""
+import fenics as fn
+
+__version__ = "0.1"
+__author__ = "Reidmen Arostica"
+__date__ = "25/06/2020"
+
+def F_(coordinates: fn.Function) -> fn.grad:
+    """
+    Computes the deformation gradient operator (dx/dX).
+    """
+    return fn.grad(coordinates)
+
+def inv_F_(coordinates: fn.Function) -> fn.inv:
+    """
+    Computes the inverse of the deformation gradient.
+    """
+    return fn.inv(F_(coordinates))
+
+def J_(coordinates: fn.Function) -> fn.det:
+    """
+    Computes the determinant of the deformation gradient.
+    """
+    return fn.det(F_(coordinates))
+
+def grad_(u: fn.Function, coordinates: fn.Function) -> fn.Function:
+    """
+    Computes the grad operator on the reference configuration.
+    """
+    inv_F = fn.inv(F_(coordinates))
+    return fn.dot(fn.grad(u), inv_F)
+
+def nabla_grad_(u: fn.Function, coordinates: fn.Function) -> fn.Function:
+    """
+    Computes the nabla_grad on the reference configuration.
+    """
+    return grad_(u, coordinates).T
+
+def sym_grad_(u: fn.Function, coordinates: fn.Function) -> fn.Function:
+    """
+    Computes the symmetric gradient operator w.r.t. the reference configuration.
+    """
+    return 0.5*(grad_(u, coordinates) + grad_(u, coordinates).T)
+
+def div_(u: fn.Function, coordinates: fn.Function) -> fn.tr:
+    """
+    Computes the Piola div operator w.r.t. the reference configuration.
+    """
+    J = J_(coordinates)
+    inv_F = inv_F_(coordinates)
+    return fn.div(J*fn.dot(inv_F, u))

main.py

@@ -0,0 +1,144 @@
+__version__ = "0.2"
+__author__ = "Reidmen Arostica"
+__date__ = "18/07/2020"
+
+# Libraries 
+from solvers import iNSE_Monolithic_Solver, iNSE_CT_Solver
+from solvers import testing_comparative_plots
+
+# Basic info to user
+# print("SOLVER VERSION: {}".format(__version__))
+# print("__author__: {}".format(__author__))
+
+
+def testing_defs_ALE(solver_to_use='LU', gcl_stab=True, jac_at_n=False):
+    """
+    Computes the iNSE problem under several deformations (including the Identity).
+    - Identity -> ('x[0]', 'x[1]')
+    - High Oscillation (init with expansion  ) -> ('x[0]*(1+0.9*sin(16*pi*t/2))', 'x[1]')
+    - High Oscillation (init with contraction) -> ('x[0]*(1+0.9*sin(16*pi*t/2))', 'x[1]')
+    """
+    print("""Initializing testing of two deformation cases with input parameters:
+    Solver to use -solver_to_use- : {}
+    GCL-stabilization term -gcl_stab- : {}
+    Evaluation of jacobian at n in time-derivate term -jac_at_n- : {}
+    """.format(solver_to_use, gcl_stab, jac_at_n))
+    # Define dict of deformation used and output names
+    dict_def = {
+    #    'Identity': ('x[0] + 0*t', 'x[1]'),
+        'High_Osc_exp': ('x[0]*(1+0.9*sin(16*pi*t/2))', 'x[1]'),
+        'High_Osc_con': ('x[0]*(1-0.9*sin(16*pi*t/2))', 'x[1]')
+    }
+    dict_def_names = {
+    #    'Identity': ('Idx', 'Idy'),
+        'High_Osc_exp': ('High_Osc_exp', 'Idy'),
+        'High_Osc_con': ('High_Osc_con', 'Idy')
+    }
+    dict_defs = {def_type: None for def_type in dict_def.keys()}
+    dict_vels = {def_type: None for def_type in dict_def.keys()}
+    dict_pres = {def_type: None for def_type in dict_def.keys()}
+    for def_type in dict_defs.keys():
+        print("\nInitializing iNSE problem with Deformation type: {}\n".format(def_type))
+        parameters = {
+            'gcl_stab': gcl_stab, # True, False whenever we want to impose CGL in the scheme
+            'jac_at_n': jac_at_n, # Choose True, False for taking the Jacobian evaluated at time k+1 or k
+            'init_cond': True,
+            'Nx': 6*16,
+            'Ny': 16,
+            'T': 2.0, #0.5, 1.0, 2.0
+            'mu': 0.035,
+            'rho_f': 1.2,
+            'tau': 0.01, # 0.01
+            'deg_U': 2,
+            'deg_P': 1
+        }
+        deformation = { 
+            'expression': dict_def[def_type],
+            'output_name': dict_def_names[def_type]
+        }
+        stabilization = {
+            'type': None
+        }
+        print("\nSetup with initial condition: {}".format(True))
+        dict_defs[def_type], dict_vels[def_type], dict_pres[def_type] = iNSE_Monolithic_Solver(parameters, deformation, stabilization, solver_to_use=solver_to_use)
+        print("\n----------------Computation Completed!----------------\n")
+    print("\n----------------Simulations Done!----------------.\n")
+
+    # Return dictionaries with values for further processing...
+    return dict_defs, dict_vels, dict_pres
+
+
+def testing_defs_ALE_CT(solver_to_use='LU', gcl_stab=True, jac_at_n=False):
+    """
+    Computes the iNSE problem under several deformations (including the Identity).
+    - Identity -> ('x[0]', 'x[1]')
+    - High Oscillation (init with expansion  ) -> ('x[0]*(1+0.9*sin(16*pi*t/2))', 'x[1]')
+    - High Oscillation (init with contraction) -> ('x[0]*(1+0.9*sin(16*pi*t/2))', 'x[1]')
+    """
+    print("""Initializing testing of two deformation cases with input parameters:
+    Solver to use -solver_to_use- : {}
+    """.format(solver_to_use))
+    # Define dict of deformation used and output names
+    dict_def = {
+    #    'Identity': ('x[0] + 0*t', 'x[1]'),
+        'High_Osc_exp': ('x[0]*(1+0.9*sin(16*pi*t/2))', 'x[1]'),
+        'High_Osc_con': ('x[0]*(1-0.9*sin(16*pi*t/2))', 'x[1]')
+    }
+    dict_def_names = {
+    #    'Identity': ('Idx', 'Idy'),
+        'High_Osc_exp': ('High_Osc_exp', 'Idy'),
+        'High_Osc_con': ('High_Osc_con', 'Idy')
+    }
+    dict_defs = {def_type: None for def_type in dict_def.keys()}
+    dict_vels = {def_type: None for def_type in dict_def.keys()}
+    dict_pres = {def_type: None for def_type in dict_def.keys()}
+    for def_type in dict_defs.keys():
+        print("\nInitializing iNSE problem with Deformation type: {}\n".format(def_type))
+        parameters = {
+            'gcl_stab': gcl_stab,
+            'jac_at_n': jac_at_n,
+            'init_cond': True,
+            'Nx': 6*16,
+            'Ny': 16,
+            'T': 2.0, #0.5, 1.0, 2.0
+            'mu': 0.035,
+            'rho_f': 1.2,
+            'tau': 0.01, # 0.01
+            'deg_U': 1,
+            'deg_P': 1
+        }
+        deformation = { 
+            'expression': dict_def[def_type],
+            'output_name': dict_def_names[def_type]
+        }
+        stabilization = {
+            'type': None
+        }
+        print("\nSetup with initial condition: {}".format(True))
+        dict_defs[def_type], dict_vels[def_type], dict_pres[def_type] = iNSE_CT_Solver(parameters, deformation, stabilization, solver_to_use=solver_to_use)
+        print("\n----------------Computation Completed!----------------\n")
+    print("\n----------------Simulations Done!----------------.\n")
+
+    # Return dictionaries with values for further processing...
+    return dict_defs, dict_vels, dict_pres
+
+
+def testing_ALE(solver_to_use='LU', gcl_stab=True):
+    """
+    Computes testing_defs_ALE (each deformation case) for two different monolithic schemes, i.e.
+    using evaluation of J at {n, n+1}, both with GCL + Temam stabilizations.
+    """
+    print("The solution of the linear system is done with built-in PETS Krylov solver.")
+    dict_fields_mt = {jac_name: None for jac_name in ['j_implicit', 'j_explicit']}
+    dict_fields_ct = {jac_name: None for jac_name in ['j_implicit', 'j_explicit']}
+    for jac_name, jac_at_n in zip(['j_implicit', 'j_explicit'], [False, True]):
+        dict_fields_mt[jac_name] = testing_defs_ALE(solver_to_use=solver_to_use, gcl_stab=gcl_stab, jac_at_n=jac_at_n)
+        dict_fields_ct[jac_name] = testing_defs_ALE_CT(solver_to_use=solver_to_use, gcl_stab=gcl_stab, jac_at_n=jac_at_n)
+    # Generate comparative figures for cases 'j_implicit', 'j_explicit'
+    testing_comparative_plots(dict_fields_mt, dict_fields_ct, gcl_stab, solver_to_use)
+    return dict_fields_mt, dict_fields_ct
+
+if __name__ == "__main__":
+    #testing_defs_ALE_CT(solver_to_use='Krylov')
+    #testing_defs_ALE(solver_to_use='LU')
+    testing_ALE(solver_to_use='LU') # 'LU', 'Krylov'