(Atomic) Boundary Condition Patches

Author: henryleberre

.github/workflows/frontier/build.sh

@@ -1,4 +1,4 @@
 #!/bin/bash
 
 . ./mfc.sh load -c f -m g
-./mfc.sh build -j 8 --gpu
+./mfc.sh test --dry-run -j 8 --gpu

.github/workflows/phoenix/test.sh

@@ -5,7 +5,7 @@ if [ "$job_device" == "gpu" ]; then
     build_opts="--gpu"
 fi
 
-./mfc.sh build -j 8 $build_opts
+./mfc.sh test --dry-run -j 8 $build_opts
 
 n_test_threads=8
 

.github/workflows/test.yml

@@ -84,7 +84,7 @@ jobs:
       - name: Build
         run:  |
           if [ '${{ matrix.intel }}' == 'true' ]; then . /opt/intel/oneapi/setvars.sh; fi
-          /bin/bash mfc.sh build -j $(nproc) --${{ matrix.debug }} --${{ matrix.mpi }} 
+          /bin/bash mfc.sh test --dry-run -j $(nproc) --${{ matrix.debug }} --${{ matrix.mpi }} 
 
       - name: Test
         run:  |

CMakeLists.txt

@@ -335,7 +335,7 @@ macro(HANDLE_SOURCES target useCommon)
 
         add_custom_command(
             OUTPUT   ${f90}
-            COMMAND  ${FYPP_EXE} -m re
+            COMMAND  ${FYPP_EXE} -m re -m itertools
                                  -I "${CMAKE_BINARY_DIR}/include/${target}"
                                  -I "${${target}_DIR}/include"
                                  -I "${common_DIR}/include"

docs/documentation/case.md

@@ -337,8 +337,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` |
+| `bc_[x,y,z]%%vel_beg[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]%%vel_end[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_n`   	         | Logical | Solving directly for the number density (in the method of classes) and compute void fraction from the number density |
@@ -463,7 +463,7 @@ The value of `dt` needs to be sufficiently small to satisfy the Courant-Friedric
 To newly start the simulation, set `t_step_start = 0`.
 To restart the simulation from $k$-th time step, set `t_step_start = k`; see [Restarting Cases](running.md#restarting-cases).
 
-##### Adaptive Time-Stepping
+##### Adaptive Time-Stepping (optional)
 
 - `cfl_adap_dt` enables adaptive time stepping with a constant CFL when true
 
@@ -480,6 +480,32 @@ To restart the simulation from $k$-th time step, set `t_step_start = k`; see [Re
 To newly start the simulation, set `n_start = 0`.
 To restart the simulation from $k$-th time step, see [Restarting Cases](running.md#restarting-cases).
 
+#### Boundary Condition Patches (optional)
+
+Boudary condition patches allow you to define boundary conditions with more granularity
+than `bc_[x,y,z]`. When using this feature, any point along the edge of the domain can
+be assigned its own boundary condition type. Since the boundaries of a 3D domain are
+2D surfaces, the concept of patches is re-introduced.
+
+Boundary conditions are applied using the [Painter's algorithm](https://en.wikipedia.org/wiki/Painter%27s_algorithm)
+where patches with higher indices take priority, in layers. The lowest priority is given
+to `bc_[x,y,z]`. This feature is opt-in and enabled by assigning `num_bc_patches` to a
+positive value.
+
+| Parameter              | Type    | Description                                    |
+| ---:                   | :----:  | :---                                           |
+| `num_bc_patches`       | Integer | Number of boundary condition patches (default 0) |
+| `%%type` *             | Integer | Boundary condition type (negative integer, see table [Boundary Conditions](#boundary-conditions)) |
+| `%%dir` *              | Integer | About the [1] x; [2] y; [3] z; axis |
+| `%%loc` *              | Integer | About the line or surface at the [-1] beginning; [+2] end; of the `%%dir` axis. |
+| `%%geometry` *         | Integer | The geometry of the boundary condition patch. See table [Boundary Condition Patch Geometry Types](#boundary-condition-patch-geometry-types) |
+| `%%vel(i)` *           | Real    | Meaning depends on the boundary condition |
+| `%%centroid(i)` *      | Real    | Meaning depends on the patch geometry. Index $i = \text{dir}$ is always ignored |
+| `%%length(i)` *        | Real    | Meaning depends on the patch geometry. Index $i = \text{dir}$ is always ignored |
+| `%%radius` *           | Real    | Meaning depends on the patch geometry. Index $i = \text{dir}$ is always ignored |
+
+*: These parameters should be prepended with `patch_bc(j)%` where $j$ is the boundary condition patch index.
+
 ### 7. Formatted Output
 
 | Parameter            | Type    | Description                                    |
@@ -862,6 +888,13 @@ This includes types exclusive to one-, two-, and three-dimensional problems.
 The patch type number (`#`) corresponds to the input value in `input.py` labeled  `patch_icpp(j)%%geometry` where $j$ is the patch index.
 Each patch requires a different set of parameters, which are also listed in this table.
 
+### Boundary Condition Patch Geometry Types
+
+| #    | Name               | Dim.   | Requirements                          |
+| ---: | :----:             | :---   | :---                                  |
+| 1    | Cuboid             | 1 & 2  | `%%centroid(1:3)` and `%%length(1:3)` |
+| 2    | Spheroid           | 1 & 2  | `%%centroid(1:3)` and `%%radius`      |
+
 ### Immersed Boundary Patch Types
 
 | #    | Name               | Dim.   |

examples/1D_impact/case.py

@@ -25,7 +25,7 @@
                     'dt'                           : mydt,                      
                     't_step_start'                 : 0,                         
                     't_step_stop'                  : int(Nt),                        
-                    't_step_save'                  : int(math.ceil(Nt/1.)),    
+                    't_step_save'                  : 1,    
 		    # ==========================================================
 
                     # Simulation Algorithm Parameters ==========================

examples/1D_inert_shocktube/export.py

@@ -0,0 +1,30 @@
+import csv
+import numpy as np
+import statistics
+from tqdm import tqdm
+
+import mfc.viz
+from case import dt, sol_L as sol
+
+
+case = mfc.viz.Case(".", dt)
+
+for name in tqdm(sol.species_names, desc="Loading Variables"):
+    case.load_variable(f"Y_{name}", f"prim.{5 + sol.species_index(name)}")
+case.load_variable("rho", "prim.1")
+case.load_variable("u", "prim.2")
+case.load_variable("p", "prim.3")
+case.load_variable("T", "prim.15")
+
+steps = case.get_timesteps()
+
+for step in [min(steps), max(steps)]:
+    t = step * dt
+
+    with open(f"mfc-{step}.csv", "w") as f:
+        writer = csv.writer(f)
+        keys = ['x'] + list(set(case.get_data()[0].keys()) - set(["x"]))
+        writer.writerow(keys)
+        for ix, x in enumerate(sorted(case.get_coords()[0])):
+            row = [case.get_data()[step][key][ix] for key in keys]
+            writer.writerow(row)

examples/1D_reactive_shocktube/export.py

@@ -0,0 +1,30 @@
+import csv
+import numpy as np
+import statistics
+from tqdm import tqdm
+
+import mfc.viz
+from case import dt, sol_L as sol
+
+
+case = mfc.viz.Case(".", dt)
+
+for name in tqdm(sol.species_names, desc="Loading Variables"):
+    case.load_variable(f"Y_{name}", f"prim.{5 + sol.species_index(name)}")
+case.load_variable("rho", "prim.1")
+case.load_variable("u", "prim.2")
+case.load_variable("p", "prim.3")
+case.load_variable("T", "prim.15")
+
+steps = case.get_timesteps()
+
+for step in [min(steps), max(steps)]:
+    t = step * dt
+
+    with open(f"mfc-{step}.csv", "w") as f:
+        writer = csv.writer(f)
+        keys = ['x'] + list(set(case.get_data()[0].keys()) - set(["x"]))
+        writer.writerow(keys)
+        for ix, x in enumerate(sorted(case.get_coords()[0])):
+            row = [case.get_data()[step][key][ix] for key in keys]
+            writer.writerow(row)

examples/2D_advection/case.py

@@ -2,6 +2,8 @@
 
 import json
 
+L = 1.0
+
 # Configuring case dictionary
 print(json.dumps({
     # Logistics ================================================
@@ -10,9 +12,9 @@
 
     # Computational Domain Parameters ==========================
     'x_domain%beg'                 : 0.E+00,
-    'x_domain%end'                 : 1.E+00,
+    'x_domain%end'                 : L,
     'y_domain%beg'                 : 0.E+00,
-    'y_domain%end'                 : 1.E+00,
+    'y_domain%end'                 : L,
     'm'                            : 99,
     'n'                            : 99,
     'p'                            : 0,
@@ -46,6 +48,21 @@
     'bc_y%end'                     : -3,
     # ==========================================================
 
+    'num_bc_patches'               : 2,
+    'patch_bc(1)%type'             : -1,
+    'patch_bc(1)%geometry'         :  1,
+    'patch_bc(1)%dir'              :  1,
+    'patch_bc(1)%loc'              : -1,
+    'patch_bc(1)%centroid(2)'      : L/2,
+    'patch_bc(1)%length(2)'        : L/2,
+
+    'patch_bc(2)%type'             : -1,
+    'patch_bc(2)%geometry'         :  1,
+    'patch_bc(2)%dir'              :  1,
+    'patch_bc(2)%loc'              : +1,
+    'patch_bc(2)%centroid(2)'      : L/2,
+    'patch_bc(2)%length(2)'        : L/2,
+
     # Formatted Database Files Structure Parameters ============
     'format'                       : 1,
     'precision'                    : 2,
@@ -60,7 +77,7 @@
     'patch_icpp(1)%length_x'       : 1.E+00,
     'patch_icpp(1)%length_y'       : 1.E+00,
     'patch_icpp(1)%vel(1)'         : 100.E+00,
-    'patch_icpp(1)%vel(2)'         : 100.E+00,
+    'patch_icpp(1)%vel(2)'         : 0.E+00,
     'patch_icpp(1)%pres'           : 1.E+05,
     'patch_icpp(1)%alpha_rho(1)'   : 1000.E+00,
     'patch_icpp(1)%alpha_rho(2)'   : 1.,

examples/2D_riemann_test/case.py

@@ -19,8 +19,8 @@
     'p'                            : 0,
     'dt'                           : 8e-05,
     't_step_start'                 : 0,
-    't_step_stop'                  : 10000,
-    't_step_save'                  : 100,
+    't_step_stop'                  : 10000//1000,
+    't_step_save'                  : 10,
     # ==========================================================
 
     # Simulation Algorithm Parameters ==========================