Velocity BCs in all coordinate directions, new example cases, nonpoly…

…tropic bc refactor (#279)

Co-authored-by: Ben Wilfong <bwilfong3@gatech.edu>
Co-authored-by: Spencer Bryngelson <sbryngelson@gmail.com>

docs/documentation/case.md

@@ -304,6 +304,8 @@ Details of implementation of viscosity in MFC can be found in [Coralic (2015)](r
 | Parameter              | Type    | Description                                    |
 | ---:                   | :----:  |          :---                                  |
 | `bc_[x,y,z]\%beg[end]` | Integer | Beginning [ending] boundary condition in the $[x,y,z]$-direction (negative integer, see table [Boundary Conditions](#boundary-conditions)) |
+| `bc_[x,y,z]\%vb[1,2,3]`‡| Real   | Velocity in the (x,1), (y, 2), (z,3) direction applied to `bc_[x,y,z]%beg` |
+| `bc_[x,y,z]\%ve[1,2,3]`‡| Real   | Velocity in the (x,1), (y, 2), (z,3) direction applied to `bc_[x,y,z]%end` |
 | `model_eqns`           | Integer | Multicomponent model: [1] $\Gamma/\Pi_\infty$; [2] 5-equation; [3] 6-equation\\%;%[4] 4-equation |
 | `alt_soundspeed` *     | Logical | Alternate sound speed and $K \nabla \cdot u$ for 5-equation model |
 | `adv_alphan`	         | Logical | Equations for all $N$ volume fractions (instead of $N-1$) |
@@ -323,6 +325,7 @@ Details of implementation of viscosity in MFC can be found in [Coralic (2015)](r
 
 - \* Options that work only with `model_eqns` $=2$.
 - † Options that work only with `cyl_coord` $=$ `False`.
+- ‡ Options that work only with `bc_[x,y,z]%[beg,end] = -15` and/or `bc_[x,y,z]%[beg,end] = -16`
 
 The table lists simulation algorithm parameters.
 The parameters are used to specify options in algorithms that are used to integrate the governing equations of the multi-component flow based on the initial condition.
@@ -332,6 +335,12 @@ Details of the simulation algorithms and implementation of the WENO scheme can b
 - `bc_[x,y,z]%[beg,end]` specifies the boundary conditions at the beginning and the end of domain boundaries in each coordinate direction by a negative integer from -1 through -12.
 See table [Boundary Conditions](#boundary-conditions) for details.
 
+- `bc_[x,y,z]\%vb[1,2,3]` specifies the velocity in the (x,1), (y,2), (z,3) direction applied to `bc_[x,y,z]%beg` when using `bc_[x,y,z]%beg = -16`.
+Tangential velocities require viscosity, `weno_avg = T`, and `bc_[x,y,z]%beg = -16` to work properly. Normal velocities require `bc_[x,y,z]\%end = -15` or `\bc_[x,y,z]\%end = -16` to work properly.
+
+- `bc_[x,y,z]\%ve[1,2,3]` specifies the velocity in the (x,1), (y,2), (z,3) direction applied to `bc_[x,y,z]%beg` when using `bc_[x,y,z]%end = -16`.
+Tangential velocities require viscosity, `weno_avg = T`, and `bc_[x,y,z]\%end = 16` to work properly. Normal velocities require `bc_[x,y,z]\%end = -15` or `\bc_[x,y,z]\%end = -16` to work properly.
+
 - `model_eqns` specifies the choice of the multi-component model that is used to formulate the dynamics of the flow using integers from 1 through 3. 
 `model_eqns` $=$ 1, 2, and 3 correspond to $\Gamma$-$\Pi_\infty$ model ([Johnsen, 2008](references.md#Johnsen08)), 5-equation model ([Allaire et al., 2002](references.md#Allaire02)), and 6-equation model ([Saurel et al., 2009](references.md#Saurel09)), respectively.
 The difference of the two models is assessed by ([Schmidmayer et al., 2019](references.md#Schmidmayer19)).

examples/2D_lid_driven_cavity/README.md

@@ -0,0 +1,16 @@
+# Lid-Driven Cavity Problem (2D)
+
+Reference: Bezgin, D. A., & Buhendwa A. B., & Adams N. A. (2022). JAX-FLUIDS: A fully-differentiable high-order computational fluid dynamics solver for compressible two-phase flows. arXiv:2203.13760
+
+Reference: Ghia, U., & Ghia, K. N., & Shin, C. T. (1982). High-re solutions for incompressible flow
+using the Navier-Stokes equations and a multigrid method. Journal of Computational Physics, 48, 387-411
+
+Video: https://youtube.com/shorts/JEP28scZrBM?feature=share
+
+## Final Condition
+
+![Final Condition](final_condition.png)
+
+## Centerline Velocities
+
+![Centerline Velocities](centerline_velocities.png)

examples/2D_lid_driven_cavity/case.py

@@ -0,0 +1,83 @@
+#!/usr/bin/env python3
+
+import json
+eps = 1e-6
+# Configuring case dictionary
+print(json.dumps({
+    # Logistics ================================================
+    'run_time_info'                : 'T',
+    # ==========================================================
+
+    # Computational Domain Parameters ==========================
+    'x_domain%beg'                 : 0.,
+    'x_domain%end'                 : 1.,
+    'y_domain%beg'                 : 0.,
+    'y_domain%end'                 : 1.,
+    'm'                            : 999,
+    'n'                            : 999,
+    'p'                            : 0,
+    'dt'                           : 1.25E-04,
+    't_step_start'                 : 0,
+    't_step_stop'                  : 1600000,
+    't_step_save'                  : 8000,
+    # ==========================================================
+
+    # Simulation Algorithm Parameters ==========================
+    'num_patches'                  : 1,
+    'model_eqns'                   : 2,
+    'alt_soundspeed'               : 'F',
+    'num_fluids'                   : 2,
+    'adv_alphan'                   : 'T',
+    'mpp_lim'                      : 'F',
+    'mixture_err'                  : 'T',
+    'time_stepper'                 : 3,
+    'weno_order'                   : 5,
+    'weno_eps'                     : 1e-16,
+    'mapped_weno'                  : 'T',
+    'weno_Re_flux'                 : 'T',
+    'mp_weno'                      : 'T',
+    'weno_avg'                     : 'T',
+    'riemann_solver'               : 2,
+    'wave_speeds'                  : 1,
+    'avg_state'                    : 2,
+    'bc_x%beg'                     : -16,
+    'bc_x%end'                     : -16,
+    'bc_y%beg'                     : -16,
+    'bc_y%end'                     : -16,
+    'bc_y%ve1'                     : 0.5,
+    # ==========================================================
+
+    # Formatted Database Files Structure Parameters ============
+    'format'                       : 1,
+    'precision'                    : 2,
+    'prim_vars_wrt'                :'T',
+    'omega_wrt(3)'                 :'T',
+    'fd_order'                     : 4,
+    'parallel_io'                  :'T',
+    # ==========================================================
+
+    # Patch 1: Base ============================================
+    'patch_icpp(1)%geometry'       : 3,
+    'patch_icpp(1)%x_centroid'     : 0.5,
+    'patch_icpp(1)%y_centroid'     : 0.5,
+    'patch_icpp(1)%length_x'       : 1.,
+    'patch_icpp(1)%length_y'       : 1.,
+    'patch_icpp(1)%vel(1)'         : 0,
+    'patch_icpp(1)%vel(2)'         : 0.,
+    'patch_icpp(1)%pres'           : 5,
+    'patch_icpp(1)%alpha_rho(1)'   : 0.5,
+    'patch_icpp(1)%alpha(1)'       : 0.5,
+    'patch_icpp(1)%alpha_rho(2)'   : 0.5,
+    'patch_icpp(1)%alpha(2)'       : 0.5,
+    # ==========================================================
+
+    # Fluids Physical Parameters ===============================
+    'fluid_pp(1)%gamma'            : 1./(1.4-1.),
+    'fluid_pp(1)%pi_inf'           : 0.,
+    'fluid_pp(1)%Re(1)'            : 1e4,
+    'fluid_pp(2)%gamma'            : 1./(1.4-1.),
+    'fluid_pp(2)%pi_inf'           : 0.,
+    'fluid_pp(2)%Re(1)'            : 1e4,
+    # ==========================================================
+}))
+# ==============================================================================

examples/2D_riemann_test/case.py

@@ -10,10 +10,10 @@
     # ==========================================================
 
     # Computational Domain Parameters ==========================
-    'x_domain%beg'                 : 0.E+00,
-    'x_domain%end'                 : 1.E+00,
-    'y_domain%beg'                 : 0.E+00,
-    'y_domain%end'                 : 1.E+00,
+    'x_domain%beg'                 : 0.,
+    'x_domain%end'                 : 1.,
+    'y_domain%beg'                 : 0.,
+    'y_domain%end'                 : 1.,
     'm'                            : 499,
     'n'                            : 499,
     'p'                            : 0,
@@ -80,8 +80,8 @@
     # ==========================================================
     # Patch 1: Base ============================================
     'patch_icpp(3)%geometry'       : 3,
-    'patch_icpp(3)%x_centroid'     : 0.9E+00,
-    'patch_icpp(3)%y_centroid'     : 0.4E+00,
+    'patch_icpp(3)%x_centroid'     : 0.9,
+    'patch_icpp(3)%y_centroid'     : 0.4,
     'patch_icpp(3)%length_x'       : 0.2,
     'patch_icpp(3)%length_y'       : 0.8,
     'patch_icpp(3)%vel(1)'         : 0,
@@ -92,8 +92,8 @@
     # ==========================================================
     # Patch 1: Base ============================================
     'patch_icpp(4)%geometry'       : 3,
-    'patch_icpp(4)%x_centroid'     : 0.9E+00,
-    'patch_icpp(4)%y_centroid'     : 0.9E+00,
+    'patch_icpp(4)%x_centroid'     : 0.9,
+    'patch_icpp(4)%y_centroid'     : 0.9,
     'patch_icpp(4)%length_x'       : 0.2,
     'patch_icpp(4)%length_y'       : 0.2,
     'patch_icpp(4)%vel(1)'         : 0,

src/common/m_derived_types.fpp

@@ -41,6 +41,12 @@ module m_derived_types
     type int_bounds_info
         integer :: beg
         integer :: end
+        real(kind(0d0)) :: vb1
+        real(kind(0d0)) :: vb2
+        real(kind(0d0)) :: vb3
+        real(kind(0d0)) :: ve1
+        real(kind(0d0)) :: ve2
+        real(kind(0d0)) :: ve3
     end type int_bounds_info
 
     !> Derived type adding beginning (beg) and end bounds info as attributes

src/common/m_helper.fpp

@@ -28,7 +28,8 @@ module m_helper
  s_swap, &
  f_cross, &
  f_create_transform_matrix, &
- f_create_bbox
+ f_create_bbox, &
+ s_print_2D_array
 
 contains
 
@@ -131,6 +132,35 @@ contains
 
     end subroutine s_comp_n_from_cons
 
+    subroutine s_print_2D_array(A, div)
+
+        real(kind(0d0)), dimension(:, :) :: A
+        integer :: i, j
+        integer :: m, n
+        real :: c
+        real, optional :: div
+
+        m = size(A, 1)
+        n = size(A, 2)
+
+        if (present(div)) then
+            c = div
+        else
+            c = 1
+        end if
+
+        print *, m, n
+
+        do i = 1, m
+            do j = 1, n
+                write (*, fmt="(F12.4)", advance="no") A(i, j)/c
+            end do
+            write (*, fmt="(A1)") " "
+        end do
+        write (*, fmt="(A1)") " "
+
+    end subroutine
+
     !> Initializes non-polydisperse bubble modeling
     subroutine s_initialize_nonpoly
         integer :: ir

src/post_process/m_checker.f90

@@ -402,6 +402,81 @@ subroutine s_check_inputs()
                              'fd_order. Exiting ...')
         end if
 
+        ! Moving Boundaries Checks: x boundaries
+        if (any((/bc_x%vb1, bc_x%vb2, bc_x%vb3/) /= 0d0)) then
+            if (bc_x%beg == 15) then
+                if (any((/bc_x%vb2, bc_x%vb3/) /= 0d0)) then
+                    call s_mpi_abort("Unsupported combination of bc_x%beg and"// &
+                                     "bc_x%vb2 or bc_x%vb3. Exiting ...")
+                end if
+            elseif (bc_x%beg /= -16) then
+                call s_mpi_abort("Unsupported combination of bc_x%beg and"// &
+                                 "bc_x%vb1, bc_x%vb2, or bc_x%vb3. Exiting...")
+            end if
+        end if
+
+        if (any((/bc_x%ve1, bc_x%ve2, bc_x%ve3/) /= 0d0)) then
+            if (bc_x%end == 15) then
+                if (any((/bc_x%ve2, bc_x%ve3/) /= 0d0)) then
+                    call s_mpi_abort("Unsupported combination of bc_x%end and"// &
+                                     "bc_x%ve2 or bc_x%ve3. Exiting ...")
+                end if
+            elseif (bc_x%end /= -16) then
+                call s_mpi_abort("Unsupported combination of bc_x%end and"// &
+                                 "bc_x%ve1, bc_x%ve2, or bc_x%ve3. Exiting...")
+            end if
+        end if
+
+        ! Moving Boundaries Checks: y boundaries
+        if (any((/bc_y%vb1, bc_y%vb2, bc_y%vb3/) /= 0d0)) then
+            if (bc_y%beg == 15) then
+                if (any((/bc_y%vb1, bc_y%vb3/) /= 0d0)) then
+                    call s_mpi_abort("Unsupported combination of bc_y%beg and"// &
+                                     "bc_y%vb1 or bc_y%vb3. Exiting ...")
+                end if
+            elseif (bc_y%beg /= -16) then
+                call s_mpi_abort("Unsupported combination of bc_y%beg and"// &
+                                 "bc_y%vb1, bc_y%vb2, or bc_y%vb3. Exiting...")
+            end if
+        end if
+
+        if (any((/bc_y%ve1, bc_y%ve2, bc_y%ve3/) /= 0d0)) then
+            if (bc_y%end == 15) then
+                if (any((/bc_y%ve1, bc_y%ve3/) /= 0d0)) then
+                    call s_mpi_abort("Unsupported combination of bc_y%end and"// &
+                                     "bc_y%ve1 or bc_y%ve3. Exiting ...")
+                end if
+            elseif (bc_y%end /= -16) then
+                call s_mpi_abort("Unsupported combination of bc_y%end and"// &
+                                 "bc_y%ve1, bc_y%ve2, or bc_y%ve3. Exiting...")
+            end if
+        end if
+
+        ! Moving Boundaries Checks: z boundaries
+        if (any((/bc_z%vb1, bc_z%vb2, bc_z%vb3/) /= 0d0)) then
+            if (bc_z%beg == 15) then
+                if (any((/bc_x%vb1, bc_x%vb2/) /= 0d0)) then
+                    call s_mpi_abort("Unsupported combination of bc_z%beg and"// &
+                                     "bc_x%vb1 or bc_x%vb1. Exiting ...")
+                end if
+            elseif (bc_z%beg /= -16) then
+                call s_mpi_abort("Unsupported combination of bc_z%beg and"// &
+                                 "bc_z%vb1, bc_z%vb2, or bc_z%vb3. Exiting...")
+            end if
+        end if
+
+        if (any((/bc_z%ve1, bc_z%ve2, bc_z%ve3/) /= 0d0)) then
+            if (bc_z%end == 15) then
+                if (any((/bc_x%ve1, bc_x%ve2/) /= 0d0)) then
+                    call s_mpi_abort("Unsupported combination of bc_z%end and"// &
+                                     "bc_z%ve2 or bc_z%ve3. Exiting ...")
+                end if
+            elseif (bc_z%end /= -16) then
+                call s_mpi_abort("Unsupported combination of bc_z%end and"// &
+                                 "bc_z%ve1, bc_z%ve2, or bc_z%ve3. Exiting...")
+            end if
+        end if
+
     end subroutine s_check_inputs
 
 end module m_checker

src/post_process/m_global_parameters.fpp

@@ -279,6 +279,13 @@ contains
         bc_y%beg = dflt_int; bc_y%end = dflt_int
         bc_z%beg = dflt_int; bc_z%end = dflt_int
 
+        #:for DIM in ['x', 'y', 'z']
+            #:for DIR in [1, 2, 3]
+                bc_${DIM}$%vb${DIR}$ = 0d0
+                bc_${DIM}$%ve${DIR}$ = 0d0
+            #:endfor
+        #:endfor
+
         ! Fluids physical parameters
         do i = 1, num_fluids_max
             fluid_pp(i)%gamma = dflt_real