Merge branch 'MFlowCode:master' into master

.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

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}

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
 
 

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

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])
 

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)

src/common/m_derived_types.f90 → 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.

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