MFlowCode · henryleberre · Oct 15, 2024
@@ -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
@@ -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
 

@@ -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:  |

@@ -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"

@@ -349,8 +349,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 |
@@ -475,7 +475,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
 
@@ -492,6 +492,36 @@ 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 and experimental)
+
+> [!WARNING]  
+> This feature is currently experimental and may not produce physicall correct
+> results when certain other features are turned on.
+
+Boundary 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                                    |
@@ -874,6 +904,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.   |

@@ -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)
@@ -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)
@@ -0,0 +1,130 @@
+#!/usr/bin/env python3
+import math
+import json
+
+# Numerical setup
+Nx      = 201  # Number of grid points in x
+Ny      = 201  # Number of grid points in y
+dx      = 1./(1.*(Nx+1))  # Grid spacing in x
+dy      = 1./(1.*(Ny+1))  # Grid spacing in y
+
+Tend    = 64E-06  # End time
+Nt      = 2000 # 2000  # Number of time steps
+mydt    = Tend/(1.*Nt)  # Time step size
+
+# Configuring case dictionary
+print(json.dumps({
+        # Logistics ================================================                 
+        'run_time_info'                : 'F',                       
+        # ==========================================================
+
+        # Computational Domain Parameters ==========================
+        'x_domain%beg'                 : 0.E+00,  # x start                    
+        'x_domain%end'                 : 2.E+00,  # x end                    
+        'y_domain%beg'                 : 0.E+00,  # y start
+        'y_domain%end'                 : 1.E+00,  # y end 
+        'm'                            : Nx,  # Number of grid points in x direction                       
+        'n'                            : Ny,  # Number of grid points in y direction                       
+        'p'                            : 0,  # Number of grid points in z (for 3D, change this)                      
+        'dt'                           : 1e-6,  # Time step size                      
+        't_step_start'                 : 0,  # Start time                         
+        't_step_stop'                  : Nt,  # End time                         
+        't_step_save'                  : 500,  # Save frequency
+        # ==========================================================
+
+        # Simulation Algorithm Parameters ==========================
+        'num_patches'                  : 1,  # Two patches                        
+        'model_eqns'                   : 2,  # Number of model equations                       
+        'alt_soundspeed'               : 'F',                      
+        'num_fluids'                   : 2,        
+        'low_Mach'                     : 0,
+        'mpp_lim'                      : 'F',                      
+        # ' mixture_err'                  : 'F',                      
+        'time_stepper'                 : 3,                                               
+        'weno_order'                   : 5,                        
+        'weno_eps'                     : 1.E-16,
+        'weno_Re_flux'                 : 'F',  
+        'weno_avg'                     : 'F',
+        'mapped_weno'                  : 'F',                     
+        'null_weights'                 : 'F',                      
+        'mp_weno'                      : 'F',                      
+        'riemann_solver'               : 1,                        
+        'wave_speeds'                  : 1,                        
+        'avg_state'                    : 2,                        
+        'bc_x%beg'                     : -3,                       
+        'bc_x%end'                     : -3,                       
+        'bc_y%beg'                     : -3,  # Boundary conditions for y direction
+        'bc_y%end'                     : -3,
+        'num_bc_patches'               : 1,
+        'patch_bc(1)%type'             : -17,
+        'patch_bc(1)%dir'              : 1,
+        'patch_bc(1)%loc'              : -1,
+        'patch_bc(1)%geometry'         : 1,
+        'patch_bc(1)%centroid(1)'      : 0, 
+        'patch_bc(1)%centroid(2)'      : 0.5,
+        'patch_bc(1)%length(2)'        : 0.26,
+        'patch_bc(1)%vel(1)'           : 10,
+        'patch_bc(1)%vel(2)'           : 0,
+        # ==========================================================
+
+        # Turning on IB ================================
+        'ib'               : 'T',
+        'num_ibs'          : 2,    
+
+        # ==========================================================
+
+        # Formatted Database Files Structure Parameters ============
+        'format'                       : 1,                        
+        'precision'                    : 2,                        
+        'prim_vars_wrt'                :'T',                       
+        'parallel_io'                  :'T',                       
+        # ==========================================================
+
+        # Patch 1 (background flow) ===================
+        'patch_icpp(1)%geometry'       : 3,  # 2D geometry                
+        'patch_icpp(1)%x_centroid'     : 1.0,  # x-center                    
+        'patch_icpp(1)%y_centroid'     : 0.5,  # y-center
+        'patch_icpp(1)%length_x'       : 2.0,   # x-length                    
+        'patch_icpp(1)%length_y'       : 1.0,   # y-length
+        'patch_icpp(1)%vel(1)'         : 0.0,
+        'patch_icpp(1)%vel(2)'         : 0.0,   # y-velocity    '100*sin(3*x*pi)'
+        'patch_icpp(1)%pres'           : 1.E5,  # Pressure                    
+        'patch_icpp(1)%alpha_rho(1)'   : 1000,  # Density    
+        'patch_icpp(1)%alpha_rho(2)'   : 0.,        
+        'patch_icpp(1)%alpha(1)'       : 1,                       
+        'patch_icpp(1)%alpha(2)'       : 0.,                          
+        'patch_icpp(1)%tau_e(1)'       : 0.0,             
+
+
+        # ==========================================================
+
+        # Patch 2 (hypo material in the center) ================
+        'patch_ib(1)%geometry'       : 3,  # 2D geometry      
+        # 'patch_ib(1)%hcid'           : 201,                    
+        'patch_ib(1)%x_centroid'     : 0.5,  # x-center                    
+        'patch_ib(1)%y_centroid'     : 0.65,  # y-center                    
+        'patch_ib(1)%length_x'       : 1.0,   # x-length                   
+        'patch_ib(1)%length_y'       : 0.04,   # y-length       
+        'patch_ib(1)%slip'           : 'T',
+
+        # ==========================================================
+
+        # Patch 3 (hypo material in the center) ================
+        'patch_ib(2)%geometry'       : 3,  # 2D geometry      
+        # 'patch_ib(1)%hcid'           : 201,                    
+        'patch_ib(2)%x_centroid'     : 0.5,  # x-center                    
+        'patch_ib(2)%y_centroid'     : 0.35,  # y-center                    
+        'patch_ib(2)%length_x'       : 1.0,   # x-length                   
+        'patch_ib(2)%length_y'       : 0.04,   # y-length       
+        'patch_ib(2)%slip'           : 'T',
+
+        # Fluids Physical Parameters ===============================
+        'fluid_pp(1)%gamma'            : 1.E+00/(6.12E+00-1.E+00),   
+        'fluid_pp(1)%pi_inf'           : 6.12E+00*3.43E+08/(6.12E+00 - 1.E+00), 
+        # 'fluid_pp(1)%G'                : 0,        
+        'fluid_pp(2)%gamma'            : 1.E+00/(1.3E+00-1.E+00),   
+        'fluid_pp(2)%pi_inf'           : 1.3E+00*2.E+08/(1.3E+00 - 1.E+00),                        
+        # 'fluid_pp(2)%G'                : 2.7E+05/(2.E+00*(1.E+00 + 0.4E+00)),     
+        'fluid_pp(2)%G'                : 1.E7,                
+        # ==========================================================
+}))