diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 72c2f3ce7..a6e4817bf 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -50,7 +50,9 @@ jobs: - name: Setup Ubuntu if: matrix.os == 'ubuntu' && matrix.intel == false - run: sudo apt install tar wget make cmake gcc g++ python3 python3-dev "openmpi-*" libopenmpi-dev + run: | + sudo apt update -y + sudo apt install -y tar wget make cmake gcc g++ python3 python3-dev "openmpi-*" libopenmpi-dev - name: Setup Ubuntu (Intel) if: matrix.os == 'ubuntu' && matrix.intel == true diff --git a/CMakeLists.txt b/CMakeLists.txt index 82cf5b1e1..5ebde586a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -290,6 +290,12 @@ function(MFC_SETUP_TARGET) add_executable(${ARGS_TARGET} ${ARGS_SOURCES}) + target_include_directories(${ARGS_TARGET} PRIVATE + "${CMAKE_SOURCE_DIR}/src/common" + "${CMAKE_SOURCE_DIR}/src/common/include" + "${CMAKE_SOURCE_DIR}/src/${ARGS_TARGET}" + "${CMAKE_SOURCE_DIR}/src/${ARGS_TARGET}/include") + string(TOUPPER "${ARGS_TARGET}" ARGS_TARGET_UPPER) target_compile_definitions( ${ARGS_TARGET} PRIVATE MFC_${CMAKE_Fortran_COMPILER_ID} diff --git a/docs/documentation/case.md b/docs/documentation/case.md index a9629da60..aec0615f3 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -139,22 +139,27 @@ the grid stretching funciton is applied and has a default value of one. ### 3. Patches -| Parameter | Type | Analytical Definition | Description | -| ---: | :----: | :----: | :--- | -| `num_patches` | Integer | Not Supported | Number of initial condition geometric patches. | -| `num_fluids` | Integer | Not Supported | Number of fluids/components present in the flow. | -| `geometry` * | Integer | Not Supported | Geometry configuration of the patch. | -| `alter_patch(i)` * | Logical | Not Supported | Alter the $i$-th patch. | -| `x[y,z]_centroid` * | Real | Not Supported | Centroid of the applied geometry in the $[x,y,z]$-direction. | -| `length_x[y,z]` * | Real | Not Supported | Length, if applicable, in the $[x,y,z]$-direction. | -| `radius` * | Real | Not Supported | Radius, if applicable, of the applied geometry. | -| `smoothen` * | Logical | Not Supported | Smoothen the applied patch. | -| `smooth_patch_id` * | Integer | Not Supported | A patch with which the applied patch is smoothened. | -| `smooth_coeff` * | Real | Not Supported | Smoothen coefficient. | -| `alpha(i)` * | Real | Supported | Volume fraction of fluid $i$. | -| `alpha_rho(i)` * | Real | Supported | Partial density of fluid $i$. | -| `pres` * | Real | Supported | Pressure. | -| `vel(i)` * | Real | Supported | Velocity in direction $i$. | +| Parameter | Type | Analytical Definition | Description | +| ---: | :----: | :----: | :--- | +| `num_patches` | Integer | Not Supported | Number of initial condition geometric patches. | +| `num_fluids` | Integer | Not Supported | Number of fluids/components present in the flow. | +| `geometry` * | Integer | Not Supported | Geometry configuration of the patch. | +| `alter_patch(i)` * | Logical | Not Supported | Alter the $i$-th patch. | +| `x[y,z]_centroid` * | Real | Not Supported | Centroid of the applied geometry in the $[x,y,z]$-direction. | +| `length_x[y,z]` * | Real | Not Supported | Length, if applicable, in the $[x,y,z]$-direction. | +| `radius` * | Real | Not Supported | Radius, if applicable, of the applied geometry. | +| `smoothen` * | Logical | Not Supported | Smoothen the applied patch. | +| `smooth_patch_id` * | Integer | Not Supported | A patch with which the applied patch is smoothened. | +| `smooth_coeff` * | Real | Not Supported | Smoothen coefficient. | +| `alpha(i)` * | Real | Supported | Volume fraction of fluid $i$. | +| `alpha_rho(i)` * | Real | Supported | Partial density of fluid $i$. | +| `pres` * | Real | Supported | Pressure. | +| `vel(i)` * | Real | Supported | Velocity in direction $i$. | +| `model%filepath` | String | Not Supported | Path to an STL or OBJ file (not all OBJs are supported). | +| `model%scale(i)` | Real | Not Supported | Model's (applied) scaling factor for component $i$. | +| `model%rotate(i)` | Real | Not Supported | Model's (applied) angle of rotation about axis $i$. | +| `model%translate(i)` | Real | Not Supported | Model's $i$-th component of (applied) translation. | +| `model%spc` | Integer | Not Supported | Number of samples per cell when discretizing the model into the grid. | *: These parameters should be prepended with `patch_icpp(j)%` where $j$ is the patch index. @@ -221,6 +226,9 @@ Optimal choice of the value of `smooth_coeff` is case-dependent and left to the - `patch_icpp(j)alpha(i)`, `patch_icpp(j)alpha_rho(i)`, `patch_icpp(j)pres`, and `texttt{patch_icpp(j)vel(i)` define for $j$-th patch the void fraction of `fluid(i)`, partial density of `fluid(i)`, the pressure, and the velocity in the $i$-th coordinate direction. These physical parameters must be consistent with fluid material's parameters defined in the next subsection. See also `adv_alphan` in table [Simulation Algorithm Parameters](#5-simulation-algorithm). +- 'model%scale', 'model%rotate` and `model%translate` define how the model should be transformed to domain-space by first +scaling by `model%scale`, then rotating about the Z, X, and Y axes (using `model%rotate`), and finally translating by `model%translate`. + ### 4. Fluid Material’s | Parameter | Type | Description | @@ -339,7 +347,7 @@ Note that `time_stepper` $=$ 3 specifies the total variation diminishing (TVD), | `schlieren_alpha(i)` | Real | Intensity of the numerical Schlieren computed via `alpha(i)` | | `probe_wrt` | Logical | Write the flow chosen probes data files for each time step | | `num_probes` | Integer | Number of probes | -| `probe(i)%[x,y,z]` | Real | Coordinates of probe $i$ | +| `probe(i)%[x,y,z]` | Real | Coordinates of probe $i$ | The table lists formatted database output parameters. The parameters define variables that are outputted from simulation and file types and formats of data as well as options for post-processing. @@ -364,17 +372,17 @@ Parallel I/O enables the use of different number of processors in each of the pr ### 7. Acoustic Source -| Parameter | Type | Description | -| ---: | :----: | :--- | -| `Monopole` | Logical | Acoustic source | -| `num_mono` | Integer | Number of acoustic sources | -| `Mono(i)%pulse` | Integer | Acoustic wave form: [1] Sine [2] Gaussian [3] Square | -| `Mono(i)%npulse` | Integer | Number of pulse cycles | -| `Mono(i)%support` | Integer | Type of the spatial support of the acoustic source : [1] 1D [2] Finite width (2D) [3] Support for finite line/patch [4] General support for 3D simulation in cartesian systems [5] Support along monopole acoustic transducer [6] Support for cylindrical coordinate system along axial-dir | -| `Mono(i)%loc(j)` | Real | $j$-th coordinate of the point that consists of $i$-th source plane | -| `Mono(i)%dir` | Real | Direction of acoustic propagation | -| `Mono(i)%mag` | Real | Pulse magnitude | -| `Mono(i)%length` | Real | Spatial pulse length | +| Parameter | Type | Description | +| ---: | :----: | :--- | +| `Monopole` | Logical | Acoustic source | +| `num_mono` | Integer | Number of acoustic sources | +| `Mono(i)%pulse` | Integer | Acoustic wave form: [1] Sine [2] Gaussian [3] Square | +| `Mono(i)%npulse` | Integer | Number of pulse cycles | +| `Mono(i)%support` | Integer | Type of the spatial support of the acoustic source : [1] 1D [2] Finite width (2D) [3] Support for finite line/patch [4] General support for 3D simulation in cartesian systems [5] Support along monopole acoustic transducer [6] Support for cylindrical coordinate system along axial-dir | +| `Mono(i)%loc(j)` | Real | $j$-th coordinate of the point that consists of $i$-th source plane | +| `Mono(i)%dir` | Real | Direction of acoustic propagation | +| `Mono(i)%mag` | Real | Pulse magnitude | +| `Mono(i)%length` | Real | Spatial pulse length | The table lists acoustic source parameters. The parameters are optionally used to define a source plane in the domain that generates an acoustic wave that propagates in a specified direction normal to the source plane (one-way acoustic source). Details of the acoustic source model can be found in [Maeda and Colonius (2017)](references.md#Maeda17). @@ -429,8 +437,6 @@ The source plane is defined in the finite region of the domain: $x\in[-\infty,\i | `sigV` | Real | Standard deviation for probability density function of bubble velocity (only when qbmm is true) | | `rhoRV` | Real | Correlation coefficient for joint probability density function of bubble radius and velocity (only when qbmm is true) | - - These options work only for gas-liquid two component flows. Component indexes are required to be 1 for liquid and 2 for gas. - \* These parameters should be pretended with patch index $1$ that is filled with liquid: `fluid_pp(1)%`. @@ -478,13 +484,13 @@ Implementation of the parameterse into the model follow [Ando (2010)](references ### 9. Velocity Field Setup -| Parameter | Type | Description | -| ---: | :----: | :--- | -| `perturb_flow` | Logical | Perturb the initlal velocity field by random noise | -| `perturb_sph` | Logical | Perturb the initial partial density by random noise | -| `perturb_sph_fluid` | Integer | Fluid component whose partial density to be perturbed | -| `vel_profile` | Logical | Set the mean streamwise velocity to hyperbolic tangent profile | -| `instability_wave` | Logical | Perturb the initial velocity field by instability waves | +| Parameter | Type | Description | +| ---: | :----: | :--- | +| `perturb_flow` | Logical | Perturb the initlal velocity field by random noise | +| `perturb_sph` | Logical | Perturb the initial partial density by random noise | +| `perturb_sph_fluid` | Integer | Fluid component whose partial density to be perturbed | +| `vel_profile` | Logical | Set the mean streamwise velocity to hyperbolic tangent profile | +| `instability_wave` | Logical | Perturb the initial velocity field by instability waves | The table lists velocity field parameters. The parameters are optionally used to define initial velocity profiles and perturbations. @@ -517,34 +523,35 @@ The table lists velocity field parameters. The parameters are optionally used to | -10 | Characteristic | Constant pressure subsonic outflow | | -11 | Characteristic | Supersonic inflow | | -12 | Characteristic | Supersonic outflow | - + The boundary condition supported by the MFC are listed in table [Boundary Conditions](#boundary-conditions). Their number (`#`) corresponds to the input value in `input.py` labeled `bc_[x,y,z]%[beg,end]` (see table [Simulation Algorithm Parameters](#5-simulation-algorithm)). The entries labeled "Characteristic." are characteristic boundary conditions based on [Thompson (1987)](references.md#Thompson87) and [Thompson (1990)](references.md#Thompson90). ### Patch types -| # | Name | Dim. | Smooth | Description | -| ---: | :----: | :---: | :---: | :--- | -| 1 | Line segment | 1 | N | Requires `x_centroid` and `x_length`. | +| # | Name | Dim. | Smooth | Description | +| ---: | :----: | :---: | :---: | :--- | +| 1 | Line segment | 1 | N | Requires `x_centroid` and `x_length`. | | 2 | Circle | 2 | Y | Requires `[x,y]_centroid` and `radius`. | -| 3 | Rectangle | 2 | N | Coordinate-aligned. Requires `[x,y]_centroid` and `[x,y]_length`. | +| 3 | Rectangle | 2 | N | Coordinate-aligned. Requires `[x,y]_centroid` and `[x,y]_length`. | | 4 | Sweep line | 2 | Y | Not coordinate aligned. Requires `[x,y]_centroid` and `normal(i)`. | | 5 | Ellipse | 2 | Y | Requires `[x,y]_centroid` and `radii(i)`. | | 6 | Vortex | 2 | N | Isentropic flow disturbance. Requires `[x,y]_centroid` and `radius`. | -| 7 | 2D analytical | 2 | N | Assigns the primitive variables as analytical functions. | +| 7 | 2D analytical | 2 | N | Assigns the primitive variables as analytical functions. | | 8 | Sphere | 3 | Y | Requires `[x,y,z]_centroid` and `radius` | | 9 | Cuboid | 3 | N | Coordinate-aligned. Requires `[x,y,z]_centroid` and `[x,y,z]_length`. | -| 10 | Cylinder | 3 | Y | Requires `[x,y,z]_centroid`, `radius`, and `[x,y,z]_length`. | -| 11 | Sweep plane | 3 | Y | Not coordinate-aligned. Requires `x[y,z]_centroid` and `normal(i)`. | -| 12 | Ellipsoid | 3 | Y | Requires `[x,y,z]_centroid` and `radii(i)`. | -| 13 | 3D analytical | 3 | N | Assigns the primitive variables as analytical functions | -| 14 | Spherical Harmonic | 3 | N | Requires `[x,y,z]_centroid`, `radius`, `epsilon`, `beta` | -| 15 | 1D analytical | 1 | N | Assigns the primitive variables as analytical functions | -| 16 | 1D bubble pulse | 1 | N | Requires `x_centroid`, `length_x` | -| 17 | Spiral | 2 | N | Requires `[x,y]_centroid` | -| 18 | 2D Varcircle | 2 | Y | Requires `[x,y]_centroid`, `radius`, and `thickness` | -| 19 | 3D Varcircle | 3 | Y | Requires `[x,y,z]_centroid`, `length_z`, `radius`, and `thickness` | +| 10 | Cylinder | 3 | Y | Requires `[x,y,z]_centroid`, `radius`, and `[x,y,z]_length`. | +| 11 | Sweep plane | 3 | Y | Not coordinate-aligned. Requires `x[y,z]_centroid` and `normal(i)`. | +| 12 | Ellipsoid | 3 | Y | Requires `[x,y,z]_centroid` and `radii(i)`. | +| 13 | 3D analytical | 3 | N | Assigns the primitive variables as analytical functions | +| 14 | Spherical Harmonic | 3 | N | Requires `[x,y,z]_centroid`, `radius`, `epsilon`, `beta` | +| 15 | 1D analytical | 1 | N | Assigns the primitive variables as analytical functions | +| 16 | 1D bubble pulse | 1 | N | Requires `x_centroid`, `length_x` | +| 17 | Spiral | 2 | N | Requires `[x,y]_centroid` | +| 18 | 2D Varcircle | 2 | Y | Requires `[x,y]_centroid`, `radius`, and `thickness` | +| 19 | 3D Varcircle | 3 | Y | Requires `[x,y,z]_centroid`, `length_z`, `radius`, and `thickness` | +| 21 | Model | 2 & 3 | Y | Imports a Model (STL/OBJ). Requires `model%filepath`. | | The patch types supported by the MFC are listed in table [Patch Types](#patch-types). This includes types exclusive to one-, two-, and three-dimensional problems. The patch type number (`#`) @@ -555,13 +562,13 @@ also listed in this table. ### Monopole supports | # | Description | -| --- | ---- | -| 1 | 1D normal to x-axis | -| 2 | 2D semi-infinite source plane | -| 3 | 3D semi-infinite source plane along some lines | -| 4 | 3D semi-infinite source plane | -| 5 | Transducer | -| 6 | Cyl_coord along axial-dir| +| ---- | ---- | +| 1 | 1D normal to x-axis | +| 2 | 2D semi-infinite source plane | +| 3 | 3D semi-infinite source plane along some lines | +| 4 | 3D semi-infinite source plane | +| 5 | Transducer | +| 6 | Cyl_coord along axial-dir | The monopole support types available in MFC are listed in table [Monopole supports](#monopole-supports). This includes types exclusive to one-, two-, and three-dimensional problems with special souce geometry like transducers as well as coordinate systems such as cylindrical coordinates. The monopole support number (`#`) corresponds to the input value in `input.py` labeled `Mono(i)%support` where @@ -569,27 +576,27 @@ $i$ is the monopole source index. ### Conservative Variables Ordering -| 5-eqn | 6-eqn | -| ---- | ---- | +| 5-eqn | 6-eqn | +| ---- | ---- | | num_fluids continuity variables | num_fluids continuity variables | -| num_dims momentum variables | num_dims momentum variables | -| 1 energy variable | 1 energy variable | -| num_fluids advection variables | num_fluids advection variables | -| N/A | num_fluids internal energy variables | +| num_dims momentum variables | num_dims momentum variables | +| 1 energy variable | 1 energy variable | +| num_fluids advection variables | num_fluids advection variables | +| N/A | num_fluids internal energy variables | The above variables are used for all simulations. -| 5-eqn | 6-eqn | -| ---- | ---- | -| sub-grid bubble variables | N/A | -| hypoelastic variables | N/A | +| 5-eqn | 6-eqn | +| ---- | ---- | +| sub-grid bubble variables | N/A | +| hypoelastic variables | N/A | The above variables correspond to optional physics. ### Primitive Variables Ordering -| 5-eqn | 6-eqn | -| ---- | ---- | +| 5-eqn | 6-eqn | +| ----- | ---- | | num_fluids densities | num_fluids densities | | num_dims velocities | num_dims velocities | | 1 pressure | 1 pressure | diff --git a/mfc.sh b/mfc.sh index 4a8d0bbc3..d778921d9 100755 --- a/mfc.sh +++ b/mfc.sh @@ -230,7 +230,7 @@ else cmake_patch=$(echo $cmake_verstr | tr '.' '\n' | sed -n 3p) cmake_version="$(printf %05d%05d%05d $cmake_major $cmake_minor $cmake_patch)" - MFC_CMAKE_MIN_VERSTR=$(cat CMakeLists.txt | head -n 1 | sed 's/[^0-9,.]*//g') + MFC_CMAKE_MIN_VERSTR=$(cat CMakeLists.txt | grep cmake_minimum_required | head -n 1 | sed 's/[^0-9,.]*//g') MFC_CMAKE_MIN_MAJOR=$(echo $MFC_CMAKE_MIN_VERSTR | tr '.' '\n' | head -n 1) MFC_CMAKE_MIN_MINOR=$(echo $MFC_CMAKE_MIN_VERSTR | tr '.' '\n' | head -n 2 | tail -n 1) MFC_CMAKE_MIN_PATCH=0 @@ -379,6 +379,8 @@ if [ ! -f "$(pwd)/build/venv/bin/activate" ]; then fi ok "Created a$MAGENTA Python$COLOR_RESET virtual environment (venv)." + + rm "$(pwd)/build/requirements.txt" > /dev/null 2>&1 || true fi diff --git a/misc/img2stl.py b/misc/img2stl.py new file mode 100644 index 000000000..88092021c --- /dev/null +++ b/misc/img2stl.py @@ -0,0 +1,84 @@ +import math, argparse + +from PIL import Image, ImageFilter + +parser = argparse.ArgumentParser(prog='img2stl', description='Turn an image into an STL file') +parser.add_argument('input', type=str, help='Path to the image file') +parser.add_argument('output', type=str, help='Path to the output OBJ file') + +VARS = vars(parser.parse_args()) + +img = Image.open(VARS['input']) + +W, H = img.size + +pixels = img.load() +for y in range(H): + for x in range(W): + if pixels[x,y][3] < 255: + pixels[x,y] = (0,0,0,255) + +img.save('temp.png') +img = img.convert('L') + +THRESHOLD = 0.1 + +with open(VARS['output'], 'w') as f: + i = 0 + pixels = img.load() + for y in range(H): + for x in range(W): + c = pixels[x,y] + + if c < 255 * THRESHOLD: continue + + area = c / 255.0 + d = (1.0 - math.sqrt(area)) / 2.0 + + def tx(_x: int, _f: bool): + if _f: _x = _x + d*(1 - _x*2) + return -1.0+2.0 * (x + _x)/float(W) + + def ty(_y: int, _f: bool): + if _f: _y = _y + d*(1 - _y*2) + return 1.0-2.0 * ((1-H/float(W))/2 + (y+_y)/float(W)) + + def tz(_z: int, _f: bool): + if _f: _z = _z + d*(1 - _z*2) + return max(1.0 / float(W), 1.0 / float(H)) * (-1.0+2.0 * _z) + + def p1(_x: int, _y: int, _z: int): + return tx(_x, True), ty(_y, True), tz(_z, False) + + for _x, _y, _z in [p1(0,0,0), p1(1,0,0), p1(0,1,0), p1(1,1,0), + p1(0,0,1), p1(1,0,1), p1(0,1,1), p1(1,1,1)]: + f.write(f'v {_x} {_y} {_z}\n') + + f.write(f'f {i+1} {i+2} {i+3}\n'); f.write(f'f {i+2} {i+3} {i+4}\n') + f.write(f'f {i+5} {i+6} {i+7}\n'); f.write(f'f {i+6} {i+7} {i+8}\n') + + i = i + 8 + + def p2(_x: int, _y: int, _z: int): + return tx(_x, False), ty(_y, True), tz(_z, True) + + for _x, _y, _z in [p2(0,0,0), p2(1,0,0), p2(0,1,0), p2(1,1,0), + p2(0,0,1), p2(1,0,1), p2(0,1,1), p2(1,1,1)]: + f.write(f'v {_x} {_y} {_z}\n') + + f.write(f'f {i+1} {i+5} {i+7}\n'); f.write(f'f {i+1} {i+3} {i+7}\n') + f.write(f'f {i+2} {i+6} {i+8}\n'); f.write(f'f {i+2} {i+4} {i+8}\n') + + i = i + 8 + + def p3(_x: int, _y: int, _z: int): + return tx(_x, True), ty(_y, False), tz(_z, True) + + for _x, _y, _z in [p3(0,0,0), p3(1,0,0), p3(0,1,0), p3(1,1,0), + p3(0,0,1), p3(1,0,1), p3(0,1,1), p3(1,1,1)]: + f.write(f'v {_x} {_y} {_z}\n') + + f.write(f'f {i+1} {i+2} {i+5}\n'); f.write(f'f {i+2} {i+5} {i+6}\n') + f.write(f'f {i+3} {i+4} {i+7}\n'); f.write(f'f {i+4} {i+7} {i+8}\n') + + i = i + 8 diff --git a/misc/viz.py b/misc/viz.py index 9cc743920..8fc32eaf0 100644 --- a/misc/viz.py +++ b/misc/viz.py @@ -114,7 +114,7 @@ def calculate_layout(n): fig, axes = plt.subplots(rows, cols, figsize=PLOT_DIMS) - for i, name in enumerate(df.keys()[1:]): + for i, name in enumerate(VARS.keys()): print(f" * Plotting {name}...") plot(name, t_step, axes[i % rows, i // rows]) diff --git a/src/common/include/macros.fpp b/src/common/include/macros.fpp index a728c2b82..1966a0b09 100644 --- a/src/common/include/macros.fpp +++ b/src/common/include/macros.fpp @@ -24,3 +24,6 @@ deallocate(${', '.join(args)}$) !$acc exit data delete(${', '.join(args)}$) #:enddef DEALLOCATE + +#define t_vec3 real(kind(0d0)), dimension(1:3) +#define t_mat4x4 real(kind(0d0)), dimension(1:4,1:4) diff --git a/src/common/m_derived_types.f90 b/src/common/m_derived_types.fpp similarity index 85% rename from src/common/m_derived_types.f90 rename to src/common/m_derived_types.fpp index 5e2c0a449..80ae7c82d 100644 --- a/src/common/m_derived_types.f90 +++ b/src/common/m_derived_types.fpp @@ -2,6 +2,8 @@ !! @file m_derived_types.f90 !! @brief Contains module m_derived_types +#:include "macros.fpp" + !> @brief This file contains the definitions of all of the custom-defined !! types used in the pre-process code. module m_derived_types @@ -59,6 +61,45 @@ module m_derived_types integer, dimension(:, :, :), allocatable :: fullmom !< Moment indices for qbmm end type bub_bounds_info + !> Defines parameters for a Model Patch + type :: ic_model_parameters + character(LEN=pathlen_max) :: filepath !< + !! Path the STL file relative to case_dir. + + t_vec3 :: translate !< + !! Translation of the STL object. + + t_vec3 :: scale !< + !! Scale factor for the STL object. + + t_vec3 :: rotate !< + !! Angle to rotate the STL object along each cartesian coordinate axis, + !! in radians. + + integer :: spc !< + !! Number of samples per cell to use when discretizing the STL object. + end type ic_model_parameters + + type :: t_triangle + real(kind(0d0)), dimension(1:3, 1:3) :: v ! Vertices of the triangle + t_vec3 :: n ! Normal vector + end type t_triangle + + type :: t_ray + t_vec3 :: o ! Origin + t_vec3 :: d ! Direction + end type t_ray + + type :: t_bbox + t_vec3 :: min ! Minimum coordinates + t_vec3 :: max ! Maximum coordinates + end type t_bbox + + type :: t_model + integer :: ntrs ! Number of triangles + type(t_triangle), allocatable :: trs(:) ! Triangles + end type t_model + !> Derived type adding initial condition (ic) patch parameters as attributes !! NOTE: The requirements for the specification of the above parameters !! are strongly dependent on both the choice of the multicomponent flow @@ -79,6 +120,8 @@ module m_derived_types !! patch geometries. It is specified through its x-, y-, and z-components !! respectively. + type(ic_model_parameters) :: model !< Model parameters + real(kind(0d0)) :: epsilon, beta !< !! The isentropic vortex parameters administrating, respectively, both !! the amplitude of the disturbance as well as its domain of influence. diff --git a/src/common/m_helper.fpp b/src/common/m_helper.fpp index 39587cdae..54424ff86 100644 --- a/src/common/m_helper.fpp +++ b/src/common/m_helper.fpp @@ -2,6 +2,7 @@ !> !! @file m_helper.f90 !! @brief Contains module m_helper + module m_helper ! Dependencies ============================================================= @@ -20,7 +21,14 @@ module m_helper s_initialize_nonpoly, & s_simpson, & s_transcoeff, & - s_int_to_str + s_int_to_str, & + s_transform_vec, & + s_transform_triangle, & + s_transform_model, & + s_swap, & + f_cross, & + f_create_transform_matrix, & + f_create_bbox contains @@ -148,15 +156,9 @@ contains rhol0 = rhoref pl0 = pref -#ifdef MFC_SIMULATION @:ALLOCATE(pb0(nb), mass_n0(nb), mass_v0(nb), Pe_T(nb)) @:ALLOCATE(k_n(nb), k_v(nb), omegaN(nb)) @:ALLOCATE(Re_trans_T(nb), Re_trans_c(nb), Im_trans_T(nb), Im_trans_c(nb)) -#else - allocate(pb0(nb), mass_n0(nb), mass_v0(nb), Pe_T(nb)) - allocate(k_n(nb), k_v(nb), omegaN(nb)) - allocate(Re_trans_T(nb), Re_trans_c(nb), Im_trans_T(nb), Im_trans_c(nb)) -#endif pb0(:) = dflt_real mass_n0(:) = dflt_real @@ -331,5 +333,150 @@ contains tmp = DEXP(-0.5d0*(phi(nb)/sd)**2)/DSQRT(2.d0*pi)/sd weight(nb) = tmp*dphi/3.d0 end subroutine s_simpson + + !> This procedure computes the cross product of two vectors. + !! @param a First vector. + !! @param b Second vector. + !! @return The cross product of the two vectors. + function f_cross(a, b) result(c) + real(kind(0d0)), dimension(3), intent(in) :: a, b + real(kind(0d0)), dimension(3) :: c + + c(1) = a(2) * b(3) - a(3) * b(2) + c(2) = a(3) * b(1) - a(1) * b(3) + c(3) = a(1) * b(2) - a(2) * b(1) + end function f_cross + + !> This procedure swaps two real numbers. + !! @param lhs Left-hand side. + !! @param rhs Right-hand side. + subroutine s_swap(lhs, rhs) + real(kind(0d0)), intent(inout) :: lhs, rhs + real(kind(0d0)) :: ltemp + + ltemp = lhs + lhs = rhs + rhs = ltemp + end subroutine s_swap + + !> This procedure creates a transformation matrix. + !! @param p Parameters for the transformation. + !! @return Transformation matrix. + function f_create_transform_matrix(p) result(out_matrix) + + type(ic_model_parameters) :: p + + t_mat4x4 :: sc, rz, rx, ry, tr, out_matrix + + sc = transpose(reshape([ & + p%scale(1), 0d0, 0d0, 0d0, & + 0d0, p%scale(2), 0d0, 0d0, & + 0d0, 0d0, p%scale(3), 0d0, & + 0d0, 0d0, 0d0, 1d0 ], shape(sc))) + + rz = transpose(reshape([ & + cos(p%rotate(3)), -sin(p%rotate(3)), 0d0, 0d0, & + sin(p%rotate(3)), cos(p%rotate(3)), 0d0, 0d0, & + 0d0, 0d0, 1d0, 0d0, & + 0d0, 0d0, 0d0, 1d0 ], shape(rz))) + + rx = transpose(reshape([ & + 1d0, 0d0, 0d0, 0d0, & + 0d0, cos(p%rotate(1)), -sin(p%rotate(1)), 0d0, & + 0d0, sin(p%rotate(1)), cos(p%rotate(1)), 0d0, & + 0d0, 0d0, 0d0, 1d0 ], shape(rx))) + + ry = transpose(reshape([ & + cos(p%rotate(2)), 0d0, sin(p%rotate(2)), 0d0, & + 0d0, 1d0, 0d0, 0d0, & + -sin(p%rotate(2)), 0d0, cos(p%rotate(2)), 0d0, & + 0d0, 0d0, 0d0, 1d0 ], shape(ry))) + + tr = transpose(reshape([ & + 1d0, 0d0, 0d0, p%translate(1), & + 0d0, 1d0, 0d0, p%translate(2), & + 0d0, 0d0, 1d0, p%translate(3), & + 0d0, 0d0, 0d0, 1d0 ], shape(tr))) + + out_matrix = matmul(tr, matmul(ry, matmul(rx, matmul(rz, sc)))) + + end function f_create_transform_matrix + + !> This procedure transforms a vector by a matrix. + !! @param vec Vector to transform. + !! @param matrix Transformation matrix. + subroutine s_transform_vec(vec, matrix) + + t_vec3, intent(inout) :: vec + t_mat4x4, intent(in) :: matrix + + real(kind(0d0)), dimension(1:4) :: tmp + + tmp = matmul(matrix, [ vec(1), vec(2), vec(3), 1d0 ]) + vec = tmp(1:3) + + end subroutine s_transform_vec + + !> This procedure transforms a triangle by a matrix, one vertex at a time. + !! @param triangle Triangle to transform. + !! @param matrix Transformation matrix. + subroutine s_transform_triangle(triangle, matrix) + + type(t_triangle), intent(inout) :: triangle + t_mat4x4, intent(in) :: matrix + + integer :: i + + real(kind(0d0)), dimension(1:4) :: tmp + + do i = 1, 3 + call s_transform_vec(triangle%v(i,:), matrix) + end do + + end subroutine s_transform_triangle + + !> This procedure transforms a model by a matrix, one triangle at a time. + !! @param model Model to transform. + !! @param matrix Transformation matrix. + subroutine s_transform_model(model, matrix) + + type(t_model), intent(inout) :: model + t_mat4x4, intent(in) :: matrix + + integer :: i + + do i = 1, size(model%trs) + call s_transform_triangle(model%trs(i), matrix) + end do + + end subroutine s_transform_model + + !> This procedure creates a bounding box for a model. + !! @param model Model to create bounding box for. + !! @return Bounding box. + function f_create_bbox(model) result(bbox) + + type(t_model), intent(in) :: model + type(t_bbox) :: bbox + + integer :: i, j + + if (size(model%trs) == 0) then + bbox%min = 0d0 + bbox%max = 0d0 + return + end if + + bbox%min = model%trs(1)%v(1,:) + bbox%max = model%trs(1)%v(1,:) + + do i = 1, size(model%trs) + do j = 1, 3 + bbox%min = min(bbox%min, model%trs(i)%v(j,:)) + bbox%max = max(bbox%max, model%trs(i)%v(j,:)) + end do + end do + + end function f_create_bbox end module m_helper diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index 3eccb108f..e8118cafb 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -132,7 +132,7 @@ contains ! Depending on model_eqns and bubbles, the appropriate procedure ! for computing pressure is targeted by the procedure pointer - if ((model_eqns /= 4) .and. (bubbles .neqv. .true.)) then + if ((model_eqns /= 4) .and. (bubbles .neqv. .true.)) then pres = (energy - dyn_p - pi_inf)/gamma else if ((model_eqns /= 4) .and. bubbles) then pres = ((energy - dyn_p)/(1.d0 - alf) - pi_inf)/gamma diff --git a/src/pre_process/m_check_patches.fpp b/src/pre_process/m_check_patches.fpp index 49fe3d59d..521e481a9 100644 --- a/src/pre_process/m_check_patches.fpp +++ b/src/pre_process/m_check_patches.fpp @@ -79,7 +79,9 @@ contains elseif (patch_icpp(i)%geometry == 19) then print *, '3d var circle' elseif (patch_icpp(i)%geometry == 20) then - call s_check_2D_TaylorGreen_vortex_patch_geometry(i) + call s_check_2D_TaylorGreen_vortex_patch_geometry(i) + elseif (patch_icpp(i)%geometry == 21) then + call s_check_model_geometry(i) else call s_mpi_abort('Unsupported choice of the '// & 'geometry of active patch '//trim(iStr)//& @@ -851,4 +853,23 @@ contains end subroutine s_check_inactive_patch_primitive_variables ! ------------ + subroutine s_check_model_geometry(patch_id) ! ------------------------------ + + integer, intent(IN) :: patch_id + + logical :: file_exists + + inquire(file=patch_icpp(patch_id)%model%filepath, exist=file_exists) + + if (.not. file_exists) then + + print '(A,I0,A)', 'Model file '//trim(patch_icpp(patch_id)%model%filepath)//& + ' requested by patch ', patch_id,' does not exist. Exiting ...' + + call s_mpi_abort() + + end if + + end subroutine s_check_model_geometry ! ----------------------------------- + end module m_check_patches diff --git a/src/pre_process/m_checker.f90 b/src/pre_process/m_checker.f90 index 842d24f55..25f2c1aaf 100644 --- a/src/pre_process/m_checker.f90 +++ b/src/pre_process/m_checker.f90 @@ -26,6 +26,14 @@ subroutine s_check_inputs() logical :: dir_check !< !! Logical variable used to test the existence of folders +#ifndef MFC_MPI + if (parallel_io .eqv. .true.) then + print '(A)', 'MFC built with --no-mpi requires parallel_io=F. ' // & + 'Exiting ...' + call s_mpi_abort() + end if +#endif + bub_fac = 0 if (bubbles .and. (num_fluids == 1)) bub_fac = 1 ! Startup checks for bubbles and bubble variables diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index ab65c78f5..bb072eccb 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/src/pre_process/m_global_parameters.fpp @@ -266,6 +266,10 @@ contains do i = 1, num_patches_max patch_icpp(i)%geometry = dflt_int + patch_icpp(i)%model%scale(:) = 1d0 + patch_icpp(i)%model%translate(:) = 0d0 + patch_icpp(i)%model%filepath(:) = ' ' + patch_icpp(i)%model%spc = 10 patch_icpp(i)%x_centroid = dflt_real patch_icpp(i)%y_centroid = dflt_real patch_icpp(i)%z_centroid = dflt_real diff --git a/src/pre_process/m_initial_condition.fpp b/src/pre_process/m_initial_condition.fpp index b06f110f6..5fab877b7 100644 --- a/src/pre_process/m_initial_condition.fpp +++ b/src/pre_process/m_initial_condition.fpp @@ -119,6 +119,10 @@ contains do i = 1, num_patches + if (proc_rank == 0) then + print*, 'Processing patch', i + end if + ! Spherical patch if (patch_icpp(i)%geometry == 8) then call s_sphere(i, patch_id_fp, q_prim_vf) @@ -151,6 +155,10 @@ contains elseif (patch_icpp(i)%geometry == 19) then call s_3dvarcircle(i, patch_id_fp, q_prim_vf) + ! 3D STL patch + elseif (patch_icpp(i)%geometry == 21) then + call s_model(i, patch_id_fp, q_prim_vf) + end if end do @@ -162,6 +170,10 @@ contains do i = 1, num_patches + if (proc_rank == 0) then + print*, 'Processing patch', i + end if + ! Circular patch if (patch_icpp(i)%geometry == 2) then call s_circle(i, patch_id_fp, q_prim_vf) @@ -196,8 +208,12 @@ contains ! TaylorGreen vortex patch elseif (patch_icpp(i)%geometry == 20) then - call s_2D_TaylorGreen_vortex(i, patch_id_fp, q_prim_vf) + call s_2D_TaylorGreen_vortex(i, patch_id_fp, q_prim_vf) + ! STL patch + elseif (patch_icpp(i)%geometry == 21) then + call s_model(i, patch_id_fp, q_prim_vf) + end if end do @@ -209,6 +225,10 @@ contains do i = 1, num_patches + if (proc_rank == 0) then + print*, 'Processing patch', i + end if + ! Line segment patch if (patch_icpp(i)%geometry == 1) then call s_line_segment(i, patch_id_fp, q_prim_vf) diff --git a/src/pre_process/m_model.fpp b/src/pre_process/m_model.fpp new file mode 100644 index 000000000..29b38c6e5 --- /dev/null +++ b/src/pre_process/m_model.fpp @@ -0,0 +1,540 @@ +!> +!! @file m_model.fpp +!! @author Henry Le Berre +!! @brief Contains module m_model + +#:include 'macros.fpp' + +module m_model + + use m_helper + use m_mpi_proxy + use m_derived_types + + use iso_c_binding, only: c_char, c_int32_t, c_int16_t, c_float + + implicit none + + private + + public :: f_model_read, s_model_write, s_model_free, f_model_is_inside + +contains + + !> This procedure reads a binary STL file. + subroutine s_read_stl_binary(filepath, model) + + character(LEN=*), intent(IN) :: filepath + type(t_model), intent(OUT) :: model + + integer :: i, j, iunit, iostat + + character(kind=c_char, len=80) :: header + integer (kind=c_int32_t) :: nTriangles + + real (kind=c_float) :: normal(3), v(3, 3) + integer(kind=c_int16_t) :: attribute + + open(newunit=iunit, file=filepath, action='READ', & + form='UNFORMATTED', status='OLD', iostat=iostat, & + access='STREAM') + + if (iostat /= 0) then + print *, "Error: could not open Binary STL file ", filepath + + call s_mpi_abort() + end if + + read(iunit, iostat=iostat) header, nTriangles + + if (iostat /= 0) then + print *, "Error: could not read header from Binary STL file ", filepath + + call s_mpi_abort() + end if + + model%ntrs = nTriangles + + allocate(model%trs(model%ntrs)) + + do i = 1, model%ntrs + read(iunit) normal(:), v(1,:), v(2,:), v(3,:), attribute + + model%trs(i)%v = v + model%trs(i)%n = normal + end do + + close(iunit) + + end subroutine s_read_stl_binary + + !> This procedure reads an ASCII STL file. + subroutine s_read_stl_ascii(filepath, model) + + character(LEN=*), intent(IN) :: filepath + type(t_model), intent(OUT) :: model + + integer :: i, j, iunit, iostat + + character(80) :: line + + open(newunit=iunit, file=filepath, action='READ', & + form='FORMATTED', status='OLD', iostat=iostat, & + access='STREAM') + + if (iostat /= 0) then + print *, "Error: could not open ASCII STL file ", filepath + + call s_mpi_abort() + end if + + model%ntrs = 0 + do + if (.not. f_read_line(iunit, line)) exit + + if (line(1:6) == "facet ") then + model%ntrs = model%ntrs + 1 + end if + end do + + allocate(model%trs(model%ntrs)) + + rewind(iunit) + + i = 1 + do + if (.not. f_read_line(iunit, line)) exit + + if (line(1:5) == "solid") cycle + if (line(1:8) == "endsolid") exit + + if (line(1:12) /= "facet normal") then + print *, "Error: expected facet normal in STL file ", filepath + + call s_mpi_abort() + end if + + call s_skip_ignored_lines(iunit) + read(line(13:), *) model%trs(i)%n + + call s_skip_ignored_lines(iunit) + read(iunit, '(A)') line + + do j = 1, 3 + if (.not. f_read_line(iunit, line)) exit + + if (line(1:6) /= "vertex") then + print *, "Error: expected vertex in STL file ", filepath + + call s_mpi_abort() + end if + + call s_skip_ignored_lines(iunit) + read(line(7:), *) model%trs(i)%v(j,:) + end do + + if (.not. f_read_line(iunit, line)) exit + if (.not. f_read_line(iunit, line)) exit + + if (line(1:8) /= "endfacet") then + print *, "Error: expected endfacet in STL file ", filepath + + call s_mpi_abort() + end if + + i = i + 1 + end do + + end subroutine s_read_stl_ascii + + !> This procedure reads an STL file. + subroutine s_read_stl(filepath, model) + + character(LEN=*), intent(IN) :: filepath + type(t_model), intent(OUT) :: model + + integer :: iunit, iostat + + character(80) :: line + + open(newunit=iunit, file=filepath, action='READ', & + form='FORMATTED', status='OLD', iostat=iostat, & + access='STREAM') + + if (iostat /= 0) then + print *, "Error: could not open STL file ", filepath + + call s_mpi_abort() + end if + + read(iunit, '(A)') line + + close(iunit) + + if (line(1:5) == "solid") then + call s_read_stl_ascii(filepath, model) + else + call s_read_stl_binary(filepath, model) + end if + + end subroutine + + !> This procedure reads an OBJ file. + subroutine s_read_obj(filepath, model) + + character(LEN=*), intent(IN) :: filepath + type(t_model), intent(OUT) :: model + + integer :: i, j, k, l, iunit, iostat, nVertices + + t_vec3, allocatable :: vertices(:,:) + + character(80) :: line + + open(newunit=iunit, file=filepath, action='READ', & + form='FORMATTED', status='OLD', iostat=iostat, & + access='STREAM') + + if (iostat /= 0) then + print *, "Error: could not open model file ", filepath + + call s_mpi_abort() + end if + + nVertices = 0 + model%ntrs = 0 + do + if (.not. f_read_line(iunit, line)) exit + + select case(line(1:2)) + case ("v ") + nVertices = nVertices + 1 + case ("f ") + model%ntrs = model%ntrs + 1 + end select + end do + + rewind(iunit) + + allocate(vertices(nVertices, 1:3)) + allocate(model%trs(model%ntrs)) + + i = 1 + j = 1 + + do + if (.not. f_read_line(iunit, line)) exit + + select case (line(1:2)) + case ("g ") + case ("vn") + case ("vt") + case ("l ") + case ("v ") + read(line(3:), *) vertices(i,:) + i = i + 1 + case ("f ") + read(line(3:), *) k, l, j + model%trs(j)%v(1,:) = vertices(k,:) + model%trs(j)%v(2,:) = vertices(l,:) + model%trs(j)%v(3,:) = vertices(j,:) + j = j + 1 + case default + print *, "Error: unknown line type in OBJ file ", filepath + print *, "Line: ", line + + call s_mpi_abort() + end select + end do + + deallocate(vertices) + + close(iunit) + + end subroutine + + !> This procedure reads a mesh from a file. + !! @param filepath Path to the file to read. + !! @return The model read from the file. + function f_model_read(filepath) result(model) + + character(LEN=*), intent(IN) :: filepath + + type(t_model) :: model + + select case (filepath(len(trim(filepath))-3:len(trim(filepath)))) + case (".stl") + call s_read_stl(filepath, model) + case (".obj") + call s_read_obj(filepath, model) + case default + print *, "Error: unknown model file format for file ", filepath + + call s_mpi_abort() + end select + + end function f_model_read + + !> This procedure writes a binary STL file. + subroutine s_write_stl(filepath, model) + + character(LEN=*), intent(IN) :: filepath + type(t_model), intent(IN) :: model + + integer :: i, j, iunit, iostat + + character(kind=c_char, len=80), parameter :: header = "Model file written by MFC." + integer (kind=c_int32_t) :: nTriangles + real (kind=c_float) :: normal(3), v(3) + integer (kind=c_int16_t) :: attribute + + open(newunit=iunit, file=filepath, action='WRITE', & + form='UNFORMATTED', iostat=iostat, access='STREAM') + + if (iostat /= 0) then + print *, "Error: could not open STL file ", filepath + + call s_mpi_abort() + end if + + nTriangles = model%ntrs + write(iunit, iostat=iostat) header, nTriangles + + if (iostat /= 0) then + print *, "Error: could not write header to STL file ", filepath + + call s_mpi_abort() + end if + + do i = 1, model%ntrs + normal = model%trs(i)%n + write(iunit) normal + + do j = 1, 3 + v = model%trs(i)%v(j,:) + write(iunit) v(:) + end do + + attribute = 0 + write(iunit) attribute + end do + + close(iunit) + + end subroutine s_write_stl + + !> This procedure writes an OBJ file. + subroutine s_write_obj(filepath, model) + + character(LEN=*), intent(IN) :: filepath + type(t_model), intent(IN) :: model + + integer :: iunit, iostat + + integer :: i, j + + open(newunit=iunit, file=filepath, action='WRITE', & + form='FORMATTED', iostat=iostat, access='STREAM') + + if (iostat /= 0) then + print *, "Error: could not open OBJ file ", filepath + + call s_mpi_abort() + end if + + write(iunit, '(A)') "# Model file written by MFC." + + do i = 1, model%ntrs + do j = 1, 3 + write(iunit, '(A, " ", (f30.20), " ", (f30.20), " ", (f30.20))') & + "v", model%trs(i)%v(j,1), model%trs(i)%v(j,2), model%trs(i)%v(j,3) + end do + + write(iunit, '(A, " ", I0, " ", I0, " ", I0)') & + "f", i*3-2, i*3-1, i*3 + end do + + close(iunit) + + end subroutine s_write_obj + + !> This procedure writes a binary STL file. + !! @param filepath Path to the file to write. + !! @param triangles Triangles to write. + subroutine s_model_write(filepath, model) + + character(LEN=*), intent(IN) :: filepath + type(t_model), intent(IN) :: model + + select case (filepath(len(trim(filepath))-3:len(trim(filepath)))) + case (".stl") + call s_write_stl(filepath, model) + case (".obj") + call s_write_obj(filepath, model) + case default + print *, "Error: unknown model file format for file ", filepath + + call s_mpi_abort() + end select + + end subroutine s_model_write + + !> This procedure frees the memory allocated for an STL mesh. + subroutine s_model_free(model) + + type(t_model), intent(INOUT) :: model + + deallocate(model%trs) + + end subroutine s_model_free + + function f_read_line(iunit, line) result(bIsLine) + + integer, intent(IN) :: iunit + character(80), intent(OUT) :: line + logical :: bIsLine + + integer :: iostat + + bIsLine = .true. + + do + read(iunit, '(A)', iostat=iostat) line + + if (iostat < 0) then + bIsLine = .false. + exit + end if + + line = adjustl(trim(line)) + + if (len(trim(line)) == 0) cycle + if (line(1:5) == "solid") cycle + if (line(1:1) == "#") cycle + + exit + end do + + end function f_read_line + + subroutine s_skip_ignored_lines(iunit) + + integer, intent(IN) :: iunit + + character(80) :: line + + if (f_read_line(iunit, line)) then + backspace(iunit) + end if + + end subroutine s_skip_ignored_lines + + !> This procedure, recursively, finds whether a point is inside an octree. + !! @param model Model to search in. + !! @param point Point to test. + !! @param spacing Space around the point to search in (grid spacing). + !! @param spc Number of samples per cell. + !! @return True if the point is inside the octree, false otherwise. + function f_model_is_inside(model, point, spacing, spc) result(fraction) + + type(t_model), intent(in) :: model + t_vec3, intent(in) :: point + t_vec3, intent(in) :: spacing + integer, intent(in) :: spc + + real(kind(0d0)) :: fraction + + type(t_ray) :: ray + integer :: i, j, nInOrOut, nHits + + real(kind(0d0)), dimension(1:spc, 1:3) :: ray_origins, ray_dirs + + do i = 1, spc + call random_number(ray_origins(i,:)) + ray_origins(i,:) = point + (ray_origins(i,:) - 0.5) * spacing(:) + + call random_number(ray_dirs(i,:)) + ray_dirs(i,:) = ray_dirs(i,:) - 0.5 + ray_dirs(i,:) = ray_dirs(i,:) / sqrt(sum(ray_dirs(i,:) * ray_dirs(i,:))) + end do + + nInOrOut = 0 + do i = 1, spc + ray%o = ray_origins(i,:) + ray%d = ray_dirs(i,:) + + nHits = 0 + do j = 1, model%ntrs + if (f_intersects_triangle(ray, model%trs(j))) then + nHits = nHits + 1 + end if + end do + + nInOrOut = nInOrOut + mod(nHits, 2) + end do + + fraction = real(nInOrOut) / real(spc) + + end function f_model_is_inside + + ! From https://www.scratchapixel.com/lessons/3d-basic-rendering/ray-tracing-rendering-a-triangle/ray-triangle-intersection-geometric-solution.html + !> This procedure checks if a ray intersects a triangle. + !! @param ray Ray. + !! @param triangle Triangle. + !! @return True if the ray intersects the triangle, false otherwise. + function f_intersects_triangle(ray, triangle) result(intersects) + + type(t_ray), intent(in) :: ray + type(t_triangle), intent(in) :: triangle + + logical :: intersects + + real(kind(0d0)) :: v0v1(3), v0v2(3), N(3), P(3), C(3), edge(3), vp(3) + real(kind(0d0)) :: area2, d, t, NdotRayDirection + + intersects = .false. + + N = triangle%n + area2 = sqrt(sum(N(:) * N(:))) + + NdotRayDirection = sum(N(:) * ray%d(:)) + + if (abs(NdotRayDirection) .lt. 0.0000001) then + return + end if + + d = - sum(N(:) * triangle%v(1,:)) + t = - (sum(N(:) * ray%o(:)) + d) / NdotRayDirection + + if (t .lt. 0) then + return + end if + + P = ray%o + t * ray%d + + edge = triangle%v(2,:) - triangle%v(1,:) + vp = P - triangle%v(1,:) + C = f_cross(edge, vp) + if (sum(N(:) * C(:)) .lt. 0) then + return + end if + + edge = triangle%v(3,:) - triangle%v(2,:) + vp = P - triangle%v(2,:) + C = f_cross(edge, vp) + if (sum(N(:) * C(:)) .lt. 0) then + return + end if + + edge = triangle%v(1,:) - triangle%v(3,:) + vp = P - triangle%v(3,:) + C = f_cross(edge, vp) + if (sum(N(:) * C(:)) .lt. 0) then + return + end if + + intersects = .true. + + end function f_intersects_triangle + +end module m_model diff --git a/src/pre_process/m_mpi_proxy.fpp b/src/pre_process/m_mpi_proxy.fpp index 3df9a7cd3..b7a31c70a 100644 --- a/src/pre_process/m_mpi_proxy.fpp +++ b/src/pre_process/m_mpi_proxy.fpp @@ -71,24 +71,26 @@ contains #:endfor do i = 1, num_patches_max - call MPI_BCAST(patch_icpp(i)%geometry, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) - call MPI_BCAST(patch_icpp(i)%smooth_patch_id, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) - - call MPI_BCAST(patch_icpp(i)%smoothen, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) - call MPI_BCAST(patch_icpp(i)%alter_patch(0), num_patches_max, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) + #:for VAR in [ 'geometry', 'smooth_patch_id', 'smoothen', & + 'alter_patch' ] + call MPI_BCAST(patch_icpp(i)%${VAR}$, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + #:endfor #:for VAR in [ 'x_centroid', 'y_centroid', 'z_centroid', & & 'length_x', 'length_y', 'length_z', 'radius', 'epsilon', & & 'beta', 'smooth_coeff', 'rho', 'p0', 'm0', 'r0', 'v0', & - & 'pres', 'gamma', 'pi_inf', ] + & 'pres', 'gamma', 'pi_inf' ] call MPI_BCAST(patch_icpp(i)%${VAR}$, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) #:endfor - call MPI_BCAST(patch_icpp(i)%normal(1), 3, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - call MPI_BCAST(patch_icpp(i)%radii(1), 3, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - call MPI_BCAST(patch_icpp(i)%vel(1), 3, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - call MPI_BCAST(patch_icpp(i)%tau_e(1), 6, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - call MPI_BCAST(patch_icpp(i)%alpha_rho(1), num_fluids_max, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - call MPI_BCAST(patch_icpp(i)%alpha(1), num_fluids_max - 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + + call MPI_BCAST(patch_icpp(i)%model%filepath, len(patch_icpp(i)%model%filepath), MPI_CHARACTER, 0, MPI_COMM_WORLD, ierr) + + #:for VAR in [ 'model%translate', 'model%scale', 'model%rotate', & + 'normal', 'radii', 'vel', 'tau_e', 'alpha_rho', 'alpha' ] + call MPI_BCAST(patch_icpp(i)%${VAR}$, size(patch_icpp(i)%${VAR}$), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + #:endfor + + call MPI_BCAST(patch_icpp(i)%model%spc, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) end do ! Fluids physical parameters diff --git a/src/pre_process/m_patches.fpp b/src/pre_process/m_patches.fpp index 1d70c168d..068e97430 100644 --- a/src/pre_process/m_patches.fpp +++ b/src/pre_process/m_patches.fpp @@ -3,12 +3,14 @@ !! @brief Contains module m_patches #:include 'case.fpp' +#:include 'macros.fpp' module m_patches ! Dependencies ============================================================= - - use m_derived_types !< Definitions of the derived types + use m_model ! Subroutine(s) related to STL files + + use m_derived_types ! Definitions of the derived types use m_global_parameters !< Definitions of the global parameters @@ -38,7 +40,8 @@ module m_patches s_sphere, & s_cuboid, & s_cylinder, & - s_sweep_plane + s_sweep_plane, & + s_model real(kind(0d0)) :: x_centroid, y_centroid, z_centroid @@ -1490,8 +1493,106 @@ contains end subroutine s_sweep_plane ! ----------------------------------------- + !> The STL patch is a 2/3D geometry that is imported from an STL file. + !! @param patch_id is the patch identifier + subroutine s_model(patch_id, patch_id_fp, q_prim_vf) ! --------------------- + + integer, intent(IN) :: patch_id + integer, intent(INOUT), dimension(0:m, 0:n, 0:p) :: patch_id_fp + type(scalar_field), dimension(1:sys_size) :: q_prim_vf + + integer :: i, j, k !< Generic loop iterators + + type(t_bbox) :: bbox + type(t_model) :: model + type(ic_model_parameters) :: params + + t_vec3 :: point + + real(kind(0d0)) :: grid_mm(1:3, 1:2) + + integer :: cell_num + integer :: ncells + + t_mat4x4 :: transform + + if (proc_rank == 0) then + print*, " * Reading model: " // trim(patch_icpp(patch_id)%model%filepath) + end if + model = f_model_read(patch_icpp(patch_id)%model%filepath) + + if (proc_rank == 0) then + print*, " * Transforming model..." + end if + transform = f_create_transform_matrix(patch_icpp(patch_id)%model) + call s_transform_model(model, transform) + + bbox = f_create_bbox(model) + + if (proc_rank == 0) then + write (*,"(A, 3(2X, F20.10))") " > Model: Min:", bbox%min(1:3) + write (*,"(A, 3(2X, F20.10))") " > Cen:", (bbox%min(1:3) + bbox%max(1:3))/2d0 + write (*,"(A, 3(2X, F20.10))") " > Max:", bbox%max(1:3) + + !call s_model_write("__out__.stl", model) + !call s_model_write("__out__.obj", model) + + grid_mm(1,:) = (/ minval(x_cc) - 0d5 * dx, maxval(x_cc) + 0d5 * dx /) + grid_mm(2,:) = (/ minval(y_cc) - 0d5 * dy, maxval(y_cc) + 0d5 * dy /) + + if (p .gt. 0) then + grid_mm(3,:) = (/ minval(z_cc) - 0d5 * dz, maxval(z_cc) + 0d5 * dz /) + else + grid_mm(3,:) = (/ 0d0, 0d0 /) + end if + + write (*,"(A, 3(2X, F20.10))") " > Domain: Min:", grid_mm(:,1) + write (*,"(A, 3(2X, F20.10))") " > Cen:", (grid_mm(:,1) + grid_mm(:,2))/2d0 + write (*,"(A, 3(2X, F20.10))") " > Max:", grid_mm(:,2) + end if + + ncells = (m+1)*(n+1)*(p+1) + do i = 0, m; do j = 0, n; do k = 0, p + + cell_num = i*(n+1)*(p+1) + j*(p+1) + (k+1) + if (proc_rank == 0 .and. mod(cell_num, ncells / 100) == 0) then + write (*,"(A, I3, A)", advance="no") & + CHAR(13) // " * Generating grid: ", & + NINT(100 * real(cell_num) / ncells), "%" + end if + + point = (/ x_cc(i), y_cc(j), 0d0 /) + if (p .gt. 0) then + point(3) = z_cc(k) + end if + + if (grid_geometry == 3) then + point = f_convert_cyl_to_cart(point) + end if + + eta = f_model_is_inside(model, point, (/ dx, dy, dz /), patch_icpp(patch_id)%model%spc) + + call s_assign_patch_primitive_variables(patch_id, i, j, k, & + eta, q_prim_vf, patch_id_fp) + + ! Note: Should probably use *eta* to compute primitive variables + ! if defining them analytically. + @:analytical() + + end do; end do; end do + + if (proc_rank == 0) then + print*, "" + print*, " * Cleaning up..." + end if + + call s_model_free(model) + + end subroutine s_model ! --------------------------------------------------- + subroutine s_convert_cylindrical_to_cartesian_coord(cyl_y, cyl_z) !$acc routine seq + real(kind(0d0)), intent(IN) :: cyl_y, cyl_z cart_y = cyl_y*sin(cyl_z) @@ -1499,8 +1600,22 @@ contains end subroutine s_convert_cylindrical_to_cartesian_coord ! -------------- + function f_convert_cyl_to_cart(cyl) result(cart) + + !$acc routine seq + + t_vec3, intent(in) :: cyl + t_vec3 :: cart + + cart = (/ cyl(1), & + cyl(2)*sin(cyl(3)), & + cyl(2)*cos(cyl(3)) /) + + end function f_convert_cyl_to_cart + subroutine s_convert_cylindrical_to_spherical_coord(cyl_x, cyl_y) !$acc routine seq + real(kind(0d0)), intent(IN) :: cyl_x, cyl_y sph_phi = atan(cyl_y/cyl_x) diff --git a/src/pre_process/p_main.f90 b/src/pre_process/p_main.f90 index 31b9f6110..9908e99a3 100644 --- a/src/pre_process/p_main.f90 +++ b/src/pre_process/p_main.f90 @@ -21,8 +21,12 @@ program p_main real(kind(0d0)) :: start, finish, time_avg, time_final real(kind(0d0)), allocatable, dimension(:) :: proc_time + call random_seed() + call s_initialize_mpi_domain() + ! Initialization of the MPI environment + call s_initialize_modules() call s_read_grid() diff --git a/toolchain/mfc.py b/toolchain/mfc.py index 4af1aefff..beb5e849e 100644 --- a/toolchain/mfc.py +++ b/toolchain/mfc.py @@ -6,17 +6,17 @@ from mfc.state import ARG from mfc.run import run from mfc.test import test +from mfc.packer import packer from mfc.common import MFC_LOGO, MFCException, quit, format_list_to_string, does_command_exist from mfc.printer import cons - def __print_greeting(): MFC_LOGO_LINES = MFC_LOGO.splitlines() max_logo_line_length = max([ len(line) for line in MFC_LOGO_LINES ]) host_line = f"{getpass.getuser()}@{platform.node()} [{platform.system()}]" targets_line = f"[bold]--targets {format_list_to_string(ARG('targets'), 'magenta', 'None')}[/bold]" - help_line = "$ ./mfc.sh \[build, run, test, clean, count] --help" + help_line = "$ ./mfc.sh \[build, run, test, clean, count, packer] --help" MFC_SIDEBAR_LINES = [ "", @@ -48,8 +48,9 @@ def __checks(): def __run(): - {"test": test.test, "run": run.run, "build": build.build, - "clean": build.clean, "bench": bench.bench, "count": count.count + {"test": test.test, "run": run.run, "build": build.build, + "clean": build.clean, "bench": bench.bench, "count": count.count, + "packer": packer.packer }[ARG("command")]() diff --git a/toolchain/mfc/args.py b/toolchain/mfc/args.py index dc97925cf..f8e145251 100644 --- a/toolchain/mfc/args.py +++ b/toolchain/mfc/args.py @@ -1,9 +1,9 @@ -import re, argparse, dataclasses +import re, os.path, argparse, dataclasses from .build import get_mfc_target_names, get_target_names, get_dependencies_names from .common import format_list_to_string from .test.test import CASES as TEST_CASES - +from .packer import packer def parse(config): from .run.engines import ENGINES @@ -19,14 +19,25 @@ def parse(config): formatter_class=argparse.ArgumentDefaultsHelpFormatter, ) - parsers = parser.add_subparsers(dest="command") + parsers = parser.add_subparsers(dest="command") + run = parsers.add_parser(name="run", help="Run a case with MFC.", formatter_class=argparse.ArgumentDefaultsHelpFormatter) + test = parsers.add_parser(name="test", help="Run MFC's test suite.", formatter_class=argparse.ArgumentDefaultsHelpFormatter) + build = parsers.add_parser(name="build", help="Build MFC and its dependencies.", formatter_class=argparse.ArgumentDefaultsHelpFormatter) + clean = parsers.add_parser(name="clean", help="Clean build artifacts.", formatter_class=argparse.ArgumentDefaultsHelpFormatter) + bench = parsers.add_parser(name="bench", help="Benchmark MFC (for CI).", formatter_class=argparse.ArgumentDefaultsHelpFormatter) + count = parsers.add_parser(name="count", help="Count LOC in MFC.", formatter_class=argparse.ArgumentDefaultsHelpFormatter) + packer = parsers.add_parser(name="packer", help="Packer utility (pack/unpack/compare)", formatter_class=argparse.ArgumentDefaultsHelpFormatter) + + packers = packer.add_subparsers(dest="packer") + pack = packers.add_parser(name="pack", help="Pack a case into a single file.", formatter_class=argparse.ArgumentDefaultsHelpFormatter) + pack.add_argument("input", metavar="INPUT", type=str, default="", help="Input file of case to pack.") + pack.add_argument("-o", "--output", metavar="OUTPUT", type=str, default=None, help="Base name of output file.") - run = parsers.add_parser(name="run", help="Run a case with MFC.", formatter_class=argparse.ArgumentDefaultsHelpFormatter) - test = parsers.add_parser(name="test", help="Run MFC's test suite.", formatter_class=argparse.ArgumentDefaultsHelpFormatter) - build = parsers.add_parser(name="build", help="Build MFC and its dependencies.", formatter_class=argparse.ArgumentDefaultsHelpFormatter) - clean = parsers.add_parser(name="clean", help="Clean build artifacts.", formatter_class=argparse.ArgumentDefaultsHelpFormatter) - bench = parsers.add_parser(name="bench", help="Benchmark MFC (for CI).", formatter_class=argparse.ArgumentDefaultsHelpFormatter) - count = parsers.add_parser(name="count", help="Count LOC in MFC.", formatter_class=argparse.ArgumentDefaultsHelpFormatter) + compare = packers.add_parser(name="compare", help="Compare two cases.", formatter_class=argparse.ArgumentDefaultsHelpFormatter) + compare.add_argument("input1", metavar="INPUT1", type=str, default=None, help="First pack file.") + compare.add_argument("input2", metavar="INPUT2", type=str, default=None, help="Second pack file.") + compare.add_argument("-rel", "--reltol", metavar="RELTOL", type=float, default=1e-12, help="Relative tolerance.") + compare.add_argument("-abs", "--abstol", metavar="ABSTOL", type=float, default=1e-12, help="Absolute tolerance.") def add_common_arguments(p, mask = None): if mask is None: @@ -117,8 +128,8 @@ def append_defaults_to_data(name: str, parser): if not key in args: args[key] = val - for a, b in [("run", run ), ("test", test ), ("build", build), - ("clean", clean), ("bench", bench), ("count", count)]: + for a, b in [("run", run ), ("test", test ), ("build", build), + ("clean", clean), ("bench", bench), ("count", count)]: append_defaults_to_data(a, b) if args["command"] is None: @@ -128,4 +139,11 @@ def append_defaults_to_data(name: str, parser): # "Slugify" the name of the job args["name"] = re.sub(r'[\W_]+', '-', args["name"]) + for e in ["input", "input1", "input2"]: + if e not in args: + continue + + if args[e] is not None: + args[e] = os.path.abspath(args[e]) + return args diff --git a/toolchain/mfc/build.py b/toolchain/mfc/build.py index 34ceb5472..6074faaa0 100644 --- a/toolchain/mfc/build.py +++ b/toolchain/mfc/build.py @@ -1,8 +1,9 @@ import os, typing, dataclasses -from .state import ARG -from .printer import cons -from . import common +from .state import ARG +from .printer import cons +from . import common +from .run.input import MFCInputFile @dataclasses.dataclass @@ -189,6 +190,9 @@ def build_target(name: str, history: typing.List[str] = None): if common.system(configure, no_exception=True) != 0: raise common.MFCException(f"Failed to configure the [bold magenta]{name}[/bold magenta] target.") + if not target.isDependency and ARG("command") == "build": + MFCInputFile("", "", {}).generate(name, bOnlyFPPs = True) + common.system(build, exception_text=f"Failed to build the [bold magenta]{name}[/bold magenta] target.") common.system(install, exception_text=f"Failed to install the [bold magenta]{name}[/bold magenta] target.") diff --git a/toolchain/mfc/packer/__init__.py b/toolchain/mfc/packer/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/toolchain/mfc/packer/errors.py b/toolchain/mfc/packer/errors.py new file mode 100644 index 000000000..63a68e98d --- /dev/null +++ b/toolchain/mfc/packer/errors.py @@ -0,0 +1,52 @@ +import dataclasses, math + +@dataclasses.dataclass(repr=False) +class Error: + absolute: float + relative: float + + def __repr__(self) -> str: + return f"abs: {self.absolute:.2E}, rel: {self.relative:.2E}" + + +def compute_error(measured: float, expected: float) -> Error: + absolute = abs(measured - expected) + + if expected != 0: + relative = absolute / abs(expected) + elif measured == expected: + relative = 0 + else: + relative = float("NaN") + + return Error(absolute, relative) + + +class AverageError: + accumulated: Error + count: int + + def __init__(self) -> None: + self.accumulated = Error(0, 0) + self.count = 0 + + def get(self) -> Error: + if self.count == 0: + return Error(0, 0) + + return Error(self.accumulated.absolute / self.count, + self.accumulated.relative / self.count) + + def push(self, error: Error) -> None: + # Do not include nans in the result + # See: compute_error() + if math.isnan(error.relative): + return + + self.accumulated.absolute += error.absolute + self.accumulated.relative += error.relative + + self.count += 1 + + def __repr__(self) -> str: + return self.get().__repr__() \ No newline at end of file diff --git a/toolchain/mfc/packer/pack.py b/toolchain/mfc/packer/pack.py new file mode 100644 index 000000000..314481613 --- /dev/null +++ b/toolchain/mfc/packer/pack.py @@ -0,0 +1,94 @@ +import dataclasses, typing, os, re, math + +from .. import common +from ..common import MFCException + +from pathlib import Path + +# This class maps to the data contained in one file in D/ +@dataclasses.dataclass(repr=False) +class PackEntry: + filepath: str + doubles: typing.List[float] + + def __repr__(self) -> str: + return f"{self.filepath} {' '.join([ str(d) for d in self.doubles ])}" + + +# This class maps to the data contained in the entirety of D/: it is tush a list +# of PackEntry classes. +class Pack: + entries: typing.Dict[str, PackEntry] + + def __init__(self): + self.entries = {} + + def __init__(self, entries: typing.List[PackEntry]): + self.entries = {} + for entry in entries: + self.set(entry) + + def find(self, filepath: str) -> PackEntry: + return self.entries.get(filepath, None) + + def set(self, entry: PackEntry): + self.entries[entry.filepath] = entry + + def save(self, filepath: str): + if filepath.endswith(".py"): + filepath = os.path.dirname(filepath) + + if os.path.isdir(filepath): + filepath = os.path.join(filepath, "pack.txt") + + if not filepath.endswith(".txt"): + filepath += ".txt" + + common.file_write(filepath, '\n'.join([ str(e) for e in sorted(self.entries.values(), key=lambda x: x.filepath) ])) + + def hash_NaNs(self) -> bool: + for entry in self.entries.values(): + for double in entry.doubles: + if math.isnan(double): + return True + + return False + + +def load(filepath: str) -> Pack: + if not os.path.isfile(filepath): + filepath = os.path.join(filepath, "pack.txt") + + entries: typing.List[PackEntry] = [] + + for line in common.file_read(filepath).splitlines(): + if common.isspace(line): + continue + + arr = line.split(' ') + + entries.append(PackEntry( + filepath=arr[0], + doubles=[ float(d) for d in arr[1:] ] + )) + + return Pack(entries) + + +def compile(casepath: str) -> typing.Tuple[Pack, str]: + entries = [] + + case_dir = os.path.dirname(casepath) if os.path.isfile(casepath) else casepath + D_dir = os.path.join(case_dir, "D") + + for filepath in list(Path(D_dir).rglob("*.dat")): + short_filepath = str(filepath).replace(f'{case_dir}', '')[1:].replace("\\", "/") + + try: + doubles = [ float(e) for e in re.sub(r"[\n\t\s]+", " ", common.file_read(filepath)).strip().split(' ') ] + except ValueError: + None, f"Failed to interpret the content of [magenta]{filepath}[/magenta] as a list of floating point numbers." + + entries.append(PackEntry(short_filepath,doubles)) + + return Pack(entries), None diff --git a/toolchain/mfc/packer/packer.py b/toolchain/mfc/packer/packer.py new file mode 100644 index 000000000..1ac1bc108 --- /dev/null +++ b/toolchain/mfc/packer/packer.py @@ -0,0 +1,63 @@ +import typing, os.path + +from ..printer import cons +from ..state import ARG, ARGS +from . import pack as _pack +from . import errors +from . import tol as packtol +from ..common import MFCException + +def load(packpath: str) -> _pack.Pack: + return _pack.load(packpath) + +def pack(casepath: str, packpath: str = None) -> typing.Tuple[_pack.Pack, str]: + if packpath is None: + packpath = casepath + + p, err = _pack.compile(casepath) + + if err is not None: + return None, err + + p.save(packpath) + + return p, None + +def compare(lhs: str = None, rhs: str = None, tol: packtol.Tolerance = None) -> typing.Tuple[errors.Error, str]: + if isinstance(lhs, str): + lhs = load(lhs) + if isinstance(rhs, str): + rhs = load(rhs) + + return packtol.compare(lhs, rhs, tol) + +def packer(): + if ARG("packer") == "pack": + if not os.path.isdir(ARG("input")): + ARGS()["input"] = os.path.dirname(ARG("input")) + + out_dir = os.path.sep.join([ARG("input"), ARG("output")]) if ARG("output") is not None else None + p, err = pack(ARG("input"), out_dir) + if err is not None: + raise MFCException(err) + elif ARG("packer") == "compare": + cons.print(f"Comparing [magenta]{os.path.relpath(ARG('input1'), os.getcwd())}[/magenta] to [magenta]{os.path.relpath(ARG('input1'), os.getcwd())}[/magenta]:") + + cons.indent() + cons.print() + + tol = packtol.Tolerance(ARG("abstol"), ARG("reltol")) + + err, msg = compare(ARG("input1"), ARG("input2"), tol) + if msg is not None: + cons.print(f"[bold red]ERROR[/bold red]: The two packs are not within tolerance ({tol}).") + cons.print(msg) + else: + cons.print(f"[bold green]OK[/bold green]: The two packs are within tolerance ({tol}).") + cons.print(f"Average error: {err}.") + + cons.print() + cons.unindent() + + else: + raise MFCException(f"Unknown packer command: {ARG('packer')}") diff --git a/toolchain/mfc/packer/tol.py b/toolchain/mfc/packer/tol.py new file mode 100644 index 000000000..68eb82aaa --- /dev/null +++ b/toolchain/mfc/packer/tol.py @@ -0,0 +1,70 @@ +import math, typing + +from ..common import MFCException + +from .pack import Pack +from .errors import compute_error, AverageError, Error + +class Tolerance(Error): + pass + + +def is_close(error: Error, tolerance: Tolerance) -> bool: + if error.absolute <= tolerance.absolute: + return True + + if math.isnan(error.relative): + return True + + if error.relative <= tolerance.relative: + return True + + return False + + +def compare(candidate: Pack, golden: Pack, tol: Tolerance) -> typing.Tuple[Error, str]: + # Keep track of the average error + avg_err = AverageError() + + # Compare entry-count + if len(candidate.entries) != len(golden.entries): + return None, "Line count does not match." + + # For every entry in the golden's pack + for gFilepath, gEntry in golden.entries.items(): + # Find the corresponding entry in the candidate's pack + cEntry = candidate.find(gFilepath) + + if cEntry == None: + return None, f"No reference to {gFilepath} in the candidate's pack." + + # Compare variable-count + if len(gEntry.doubles) != len(cEntry.doubles): + return None, f"Variable count didn't match for {gFilepath}." + + # Check if each variable is within tolerance + for valIndex, (gVal, cVal) in enumerate(zip(gEntry.doubles, cEntry.doubles)): + # Keep track of the error and average errors + error = compute_error(cVal, gVal) + avg_err.push(error) + + def raise_err(msg: str): + return None, f"""\ +Variable n°{valIndex+1} (1-indexed) in {gFilepath} {msg}: + - Candidate: {cVal} + - Golden: {gVal} + - Error: {error} + - Tolerance: {tol} +""" + + if math.isnan(gVal): + raise_err("is NaN in the golden file") + + if math.isnan(cVal): + raise_err("is NaN in the pack file") + + if not is_close(error, tol): + raise_err("is not within tolerance") + + # Return the average relative error + return avg_err.get(), None diff --git a/toolchain/mfc/run/case_dicts.py b/toolchain/mfc/run/case_dicts.py index a9b09d5a3..49b37dfbf 100644 --- a/toolchain/mfc/run/case_dicts.py +++ b/toolchain/mfc/run/case_dicts.py @@ -42,6 +42,14 @@ "pi_inf", "r0", "v0", "p0", "m0"]: PRE_PROCESS.append(f"patch_icpp({p_id})%{attribute}") + PRE_PROCESS.append(f"patch_icpp({p_id})%model%filepath") + + for attribute in ["translate", "scale", "rotate"]: + for j in range(1, 4): + PRE_PROCESS.append(f"patch_icpp({p_id})%model%{attribute}({j})") + + PRE_PROCESS.append(f"patch_icpp({p_id})%model%spc") + for cmp_id, cmp in enumerate(["x", "y", "z"]): cmp_id += 1 PRE_PROCESS.append(f'patch_icpp({p_id})%{cmp}_centroid') diff --git a/toolchain/mfc/run/input.py b/toolchain/mfc/run/input.py index 1816e9e91..f38b74c47 100644 --- a/toolchain/mfc/run/input.py +++ b/toolchain/mfc/run/input.py @@ -16,6 +16,9 @@ class MFCInputFile: case_dirpath: str case_dict: dict + def __get_ndims(self) -> int: + return 1 + min(int(self.case_dict.get("n", 0)), 1) + min(int(self.case_dict.get("p", 0)), 1) + def __is_ic_analytical(self, key: str, val: str) -> bool: if common.is_number(val) or not isinstance(val, str): return False @@ -47,7 +50,6 @@ def __generate_inp(self, target_name: str) -> None: dict_str += f"{key} = {val}\n" else: dict_str += f"{key} = '{val}'\n" - continue else: ignored.append(key) @@ -86,21 +88,12 @@ def __generate_pre_fpp(self) -> None: cons.print(f"Generating [magenta]pre_process/include/case.fpp[/magenta].") cons.indent() - DIMS = { - 1: {'ptypes': [1, 15, 16], 'sf_idx': 'i, 0, 0'}, - 2: {'ptypes': [2, 3, 4, 5, 6, 7, 17, 18], 'sf_idx': 'i, j, 0'}, - 3: {'ptypes': [8, 9, 10, 11, 12, 13, 14, 19], 'sf_idx': 'i, j, k'} - } - - def get_patchtype_ndims(ptype: int) -> int: - for ndim, info in DIMS.items(): - if ptype in info['ptypes']: - return ndim - - raise common.MFCException(f"Patch of type {ptype} cannot be analytically defined.") - - PTYPES = DIMS[1]['ptypes'] + DIMS[2]['ptypes'] + DIMS[3]['ptypes'] - + DATA = { + 1: {'ptypes': [1, 15, 16], 'sf_idx': 'i, 0, 0'}, + 2: {'ptypes': [2, 3, 4, 5, 6, 7, 17, 18, 21], 'sf_idx': 'i, j, 0'}, + 3: {'ptypes': [8, 9, 10, 11, 12, 13, 14, 19, 21], 'sf_idx': 'i, j, k'} + }[self.__get_ndims()] + patches = {} for key, val in self.case_dict.items(): @@ -119,7 +112,7 @@ def get_patchtype_ndims(ptype: int) -> int: for pid, items in patches.items(): ptype = self.case_dict[f"patch_icpp({pid})%geometry"] - if ptype not in PTYPES: + if ptype not in DATA['ptypes']: raise common.MFCException(f"Patch #{pid} of type {ptype} cannot be analytically defined.") def rhs_replace(match): @@ -146,10 +139,9 @@ def rhs_replace(match): if idx != 0: qpvf_idx_offset = f" + {idx}" - ndims = get_patchtype_ndims(ptype) - sf_idx = DIMS[ndims]['sf_idx'] + sf_idx = DATA['sf_idx'] - cons.print(f"[yellow]INFO:[/yellow] {ndims}D Analytical Patch #{pid}: Code generation for [magenta]{varname}[/magenta]...") + cons.print(f"[yellow]INFO:[/yellow] {self.__get_ndims()}D Analytical Patch #{pid}: Code generation for [magenta]{varname}[/magenta]...") lhs = f"q_prim_vf({qpvf_idx_var}{qpvf_idx_offset})%sf({sf_idx})" rhs = re.sub(r"[a-zA-Z]+", rhs_replace, expr) @@ -221,8 +213,9 @@ def __generate_post_fpp(self) -> None: pass # Generate case.fpp & [target_name].inp - def generate(self, target_name: str) -> None: - self.__generate_inp(target_name) + def generate(self, target_name: str, bOnlyFPPs = False) -> None: + if not bOnlyFPPs: + self.__generate_inp(target_name) cons.print() diff --git a/toolchain/mfc/test/pack.py b/toolchain/mfc/test/pack.py deleted file mode 100644 index b46546c78..000000000 --- a/toolchain/mfc/test/pack.py +++ /dev/null @@ -1,199 +0,0 @@ -import os, re, math, typing, dataclasses - -from pathlib import Path - -from . import case -from .. import common - -from ..common import MFCException - - -@dataclasses.dataclass(repr=False) -class Error: - absolute: float - relative: float - - def __repr__(self) -> str: - return f"abs: {self.absolute:.2E}, rel: {self.relative:.2E}" - - -def compute_error(measured: float, expected: float) -> Error: - absolute = abs(measured - expected) - - if expected != 0: - relative = absolute / abs(expected) - elif measured == expected: - relative = 0 - else: - relative = float("NaN") - - return Error(absolute, relative) - - -class AverageError: - accumulated: Error - count: int - - def __init__(self) -> None: - self.accumulated = Error(0, 0) - self.count = 0 - - def get(self) -> Error: - if self.count == 0: - return Error(0, 0) - - return Error(self.accumulated.absolute / self.count, - self.accumulated.relative / self.count) - - def push(self, error: Error) -> None: - # Do not include nans in the result - # See: compute_error() - if math.isnan(error.relative): - return - - self.accumulated.absolute += error.absolute - self.accumulated.relative += error.relative - - self.count += 1 - - def __repr__(self) -> str: - return self.get().__repr__() - - -# This class maps to the data contained in one file in D/ -@dataclasses.dataclass(repr=False) -class PackEntry: - filepath: str - doubles: typing.List[float] - - def __repr__(self) -> str: - return f"{self.filepath} {' '.join([ str(d) for d in self.doubles ])}" - - -# This class maps to the data contained in the entirety of D/: it is tush a list -# of PackEntry classes. -class Pack: - entries: typing.Dict[str, PackEntry] - - def __init__(self): - self.entries = {} - - def __init__(self, entries: typing.List[PackEntry]): - self.entries = {} - for entry in entries: - self.set(entry) - - def find(self, filepath: str) -> PackEntry: - return self.entries.get(filepath, None) - - def set(self, entry: PackEntry): - self.entries[entry.filepath] = entry - - def save(self, filepath: str): - common.file_write(filepath, '\n'.join([ str(e) for e in sorted(self.entries.values(), key=lambda x: x.filepath) ])) - - -def load(filepath: str) -> Pack: - entries: typing.List[PackEntry] = [] - - for line in common.file_read(filepath).splitlines(): - if common.isspace(line): - continue - - arr = line.split(' ') - - entries.append(PackEntry( - filepath=arr[0], - doubles=[ float(d) for d in arr[1:] ] - )) - - return Pack(entries) - - -def generate(case: case.TestCase) -> Pack: - entries = [] - - case_dir = case.get_dirpath() - D_dir = os.path.join(case_dir, "D") - - for filepath in list(Path(D_dir).rglob("*.dat")): - short_filepath = str(filepath).replace(f'{case_dir}', '')[1:].replace("\\", "/") - - try: - doubles = [ float(e) for e in re.sub(r"[\n\t\s]+", " ", common.file_read(filepath)).strip().split(' ') ] - except ValueError: - raise MFCException(f"Test {case}: Failed to interpret the content of [magenta]{filepath}[/magenta] as a list of floating point numbers.") - - for double in doubles: - if math.isnan(double): - raise MFCException(f"Test {case}: A NaN was found in {filepath} while generating a pack file.") - - entries.append(PackEntry(short_filepath,doubles)) - - return Pack(entries) - - -class Tolerance(Error): - pass - - -def is_close(error: Error, tolerance: Tolerance) -> bool: - if error.absolute <= tolerance.absolute: - return True - - if math.isnan(error.relative): - return True - - if error.relative <= tolerance.relative: - return True - - return False - - -def check_tolerance(case: case.TestCase, candidate: Pack, golden: Pack, tol: float) -> Error: - # Keep track of the average error - avg_err = AverageError() - - # Compare entry-count - if len(candidate.entries) != len(golden.entries): - raise MFCException(f"Test {case}: Line count does not match.") - - # For every entry in the golden's pack - for gFilepath, gEntry in golden.entries.items(): - # Find the corresponding entry in the candidate's pack - cEntry = candidate.find(gFilepath) - - if cEntry == None: - raise MFCException(f"Test {case}: No reference of {gFilepath} in the candidate's pack.") - - # Compare variable-count - if len(gEntry.doubles) != len(cEntry.doubles): - raise MFCException(f"Test {case}: Variable count didn't match for {gFilepath}.") - - # Check if each variable is within tolerance - for valIndex, (gVal, cVal) in enumerate(zip(gEntry.doubles, cEntry.doubles)): - # Keep track of the error and average errors - error = compute_error(cVal, gVal) - avg_err.push(error) - - def raise_err(msg: str): - raise MFCException(f"""\ -Test {case}: Variable n°{valIndex+1} (1-indexed) in {gFilepath} {msg}: - - Description: {case.trace} - - Candidate: {cVal} - - Golden: {gVal} - - Error: {error} - - Tolerance: {tol} -""") - - if math.isnan(gVal): - raise_err("is NaN in the golden file") - - if math.isnan(cVal): - raise_err("is NaN in the pack file") - - if not is_close(error, Tolerance(absolute=tol, relative=tol)): - raise_err("is not within tolerance") - - # Return the average relative error - return avg_err.get() diff --git a/toolchain/mfc/test/test.py b/toolchain/mfc/test/test.py index e16f499ff..159a355fc 100644 --- a/toolchain/mfc/test/test.py +++ b/toolchain/mfc/test/test.py @@ -8,7 +8,9 @@ from .. import sched from ..common import MFCException, does_command_exist, format_list_to_string, get_program_output from ..build import build_targets, get_install_dirpath -from . import pack as packer + +from ..packer import tol as packtol +from ..packer import packer import rich, rich.table @@ -144,10 +146,14 @@ def handle_case(test: TestCase): if cmd.returncode != 0: cons.print(cmd.stdout) - raise MFCException(f"""Test {test}: Failed to execute MFC. You can find the run's output in {out_filepath}, and the case dictionary in {os.path.join(test.get_dirpath(), "case.py")}.""") + raise MFCException(f"Test {test}: Failed to execute MFC.") + + pack, err = packer.pack(test.get_dirpath()) + if err is not None: + raise MFCException(f"Test {test}: {err}") - pack = packer.generate(test) - pack.save(os.path.join(test.get_dirpath(), "pack.txt")) + if pack.hash_NaNs(): + raise MFCException(f"Test {test}: NaNs detected in the case.") golden_filepath = os.path.join(test.get_dirpath(), "golden.txt") if ARG("generate"): @@ -166,7 +172,9 @@ def handle_case(test: TestCase): golden.save(golden_filepath) else: - packer.check_tolerance(test, pack, packer.load(golden_filepath), tol) + err, msg = packtol.compare(pack, packer.load(golden_filepath), packtol.Tolerance(tol, tol)) + if msg is not None: + raise MFCException(f"Test {test}: {msg}") if ARG("test_all"): test.delete_output()