Skip to content

Commit

Permalink
(Atomic) Boundary Condition Patches
Browse files Browse the repository at this point in the history
  • Loading branch information
henryleberre committed Nov 26, 2024
1 parent e362dc9 commit c393a60
Show file tree
Hide file tree
Showing 40 changed files with 2,115 additions and 6,319 deletions.
2 changes: 1 addition & 1 deletion .github/workflows/frontier/build.sh
Original file line number Diff line number Diff line change
@@ -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
2 changes: 1 addition & 1 deletion .github/workflows/phoenix/test.sh
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
2 changes: 1 addition & 1 deletion .github/workflows/test.yml
Original file line number Diff line number Diff line change
Expand Up @@ -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: |
Expand Down
2 changes: 1 addition & 1 deletion CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down
43 changes: 40 additions & 3 deletions docs/documentation/case.md
Original file line number Diff line number Diff line change
Expand Up @@ -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 |
Expand Down Expand Up @@ -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

Expand All @@ -480,6 +480,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 |
Expand Down Expand Up @@ -862,6 +892,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. |
Expand Down
30 changes: 30 additions & 0 deletions examples/1D_inert_shocktube/export.py
Original file line number Diff line number Diff line change
@@ -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)
30 changes: 30 additions & 0 deletions examples/1D_reactive_shocktube/export.py
Original file line number Diff line number Diff line change
@@ -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)
130 changes: 130 additions & 0 deletions examples/hypotests/case.py
Original file line number Diff line number Diff line change
@@ -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,
# ==========================================================
}))
Loading

0 comments on commit c393a60

Please sign in to comment.