Upstream CheMFC infrastructure

.gitattributes

@@ -1 +1,2 @@
 /tests/**/* linguist-generated=true
+/toolchain/mechanisms/* linguist-generated=true

.github/workflows/docs.yml

@@ -1,18 +1,12 @@
 name: Documentation
 
-on:
-  push:
-    branches:
-      - master
-
-  workflow_dispatch:
+on: [push, pull_request, workflow_dispatch]
 
 jobs:
   docs:
     name:    Build & Publish
     runs-on: ubuntu-latest
 
-    if: github.repository == 'MFlowCode/MFC'
     concurrency:
       group: docs-publish
       cancel-in-progress: true
@@ -53,6 +47,7 @@ jobs:
         echo "excluded-count = ${{ steps.sitemap.outputs.excluded-count }}"
 
     - name: Publish Documentation
+      if:   github.repository == 'MFlowCode/MFC' && github.ref == 'refs/heads/master' && github.event_name == 'push'
       run:  |
         set +e
         git ls-remote "${{ secrets.DOC_PUSH_URL }}" -q

.github/workflows/test.yml

@@ -49,8 +49,9 @@ jobs:
       - name: Setup MacOS
         if:   matrix.os == 'macos'
         run:  |
-          brew install coreutils python cmake fftw hdf5 gcc@14 open-mpi
+          brew install coreutils python cmake fftw hdf5 gcc@14 boost open-mpi
           echo "FC=gfortran-14" >> $GITHUB_ENV
+          echo "BOOST_INCLUDE=/opt/homebrew/include/" >> $GITHUB_ENV
 
       - name: Setup Ubuntu
         if:   matrix.os == 'ubuntu' && matrix.intel == false
@@ -138,4 +139,3 @@ jobs:
         with:
           name: logs
           path: test-${{ matrix.device }}.out
-

CMakeLists.txt

@@ -340,6 +340,7 @@ macro(HANDLE_SOURCES target useCommon)
                                  -D MFC_${${target}_UPPER}
                                  -D MFC_COMPILER="${CMAKE_Fortran_COMPILER_ID}"
                                  -D MFC_CASE_OPTIMIZATION=False
+                                 -D chemistry=False
                                  --line-numbering
                                  --no-folding
                                  "${fpp}" "${f90}"

README.md

@@ -58,14 +58,20 @@ It's rather straightforward.
 We'll give a brief intro. here for MacOS.
 Using [brew](https://brew.sh), install MFC's dependencies:
 ```shell
-brew install coreutils python cmake fftw hdf5 gcc open-mpi
+brew install coreutils python cmake fftw hdf5 gcc boost open-mpi
 ```
 You're now ready to build and test MFC!
 Put it to a convenient directory via
 ```shell
 git clone https://github.com/MFlowCode/MFC
 cd MFC
 ```
+and be sure MFC knows where to find Boost by appending to your dotfiles and sourcing them again
+```shell
+echo -e 'export BOOST_INCLUDE=/opt/homebrew/' | tee -a ~/.bash_profile ~/.zshrc
+. ~/.bash_profile 2>/dev/null || . ~/.zshrc 2>/dev/null
+! [ -z "${BOOST_INCLUDE+x}" ] && echo 'Environment is ready!' || echo 'Error: $BOOST_INCLUDE is unset. Please adjust the previous commands to fit with your environment.'
+```
 then you can build MFC and run the test suite!
 ```shell
 ./mfc.sh build -j $(nproc)

docs/documentation/getting-started.md

@@ -101,10 +101,14 @@ You will also have access to the `.sln` Microsoft Visual Studio solution files f
 <details>
   <summary><h3>MacOS</h3></summary>
 
-Using [Homebrew](https://brew.sh/) you can install the necessary dependencies:
+Using [Homebrew](https://brew.sh/) you can install the necessary dependencies
+before configuring your environment:
 
 ```shell
-brew install coreutils python cmake fftw hdf5 gcc open-mpi
+brew install coreutils python cmake fftw hdf5 gcc boost open-mpi
+echo -e 'export BOOST_INCLUDE=/opt/homebrew/' | tee -a ~/.bash_profile ~/.zshrc
+. ~/.bash_profile 2>/dev/null || . ~/.zshrc 2>/dev/null
+! [ -z "${BOOST_INCLUDE+x}" ] && echo 'Environment is ready!' || echo 'Error: $BOOST_INCLUDE is unset. Please adjust the previous commands to fit with your environment.'
 ```
 
 They will download the dependencies MFC requires to build itself.

misc/fpp_to_fypp.sh

@@ -0,0 +1,6 @@
+#!/bin/bash
+
+for file in $(find src -type f | grep -Ev 'autogen' | grep -E '\.fpp$'); do
+    echo "$file"
+    mv "$file" "$(echo "$file" | sed s/\.fpp/\.fypp/)"
+done

src/common/include/macros.fpp

@@ -140,17 +140,19 @@
 #endif
 #:enddef
 
-#:def PROHIBIT(*args)
-    #:set condition = args[0]
-    #:if len(args) == 1
-        #:set message = '""'
-    #:else
-        #:set message = args[1]
-    #:endif
+#:def PROHIBIT(condition, message = None)
     if (${condition}$) then
-        call s_prohibit_abort("${condition}$", ${message}$)
+        call s_prohibit_abort("${condition}$", ${message or '""'}$)
     end if
 #:enddef
 
 #define t_vec3   real(kind(0d0)), dimension(1:3)
 #define t_mat4x4 real(kind(0d0)), dimension(1:4,1:4)
+
+#:def ASSERT(predicate, message = None)
+    if (.not. (${predicate}$)) then
+        call s_mpi_abort("${_FILE_.split('/')[-1]}$:${_LINE_}$: "// &
+                         "Assertion failed: ${predicate}$. " &
+                         //${message or '"No error description."'}$)
+    end if
+#:enddef

src/common/m_derived_types.fpp

@@ -8,7 +8,8 @@
 !!              types used in the pre-process code.
 module m_derived_types
 
-    use m_constants !< Constants
+    use m_constants  !< Constants
+    use m_thermochem !< Thermodynamic properties
 
     implicit none
 
@@ -192,6 +193,7 @@ module m_derived_types
         !! id for hard coded initial condition
 
         real(kind(0d0)) :: cf_val !! color function value
+        real(kind(0d0)) :: Y(1:num_species)
 
     end type ic_patch_parameters
 
@@ -300,4 +302,18 @@ module m_derived_types
 
     end type ghost_point
 
+    !> Species parameters
+    type species_parameters
+        character(LEN=name_len) :: name !< Name of species
+    end type species_parameters
+
+    !> Chemistry parameters
+    type chemistry_parameters
+        character(LEN=name_len) :: cantera_file !< Path to Cantera file
+
+        logical :: advection
+        logical :: diffusion
+        logical :: reactions
+    end type chemistry_parameters
+
 end module m_derived_types

src/common/m_finite_differences.fpp

@@ -0,0 +1,127 @@
+module m_finite_differences
+
+    use m_global_parameters
+
+    implicit none
+
+contains
+
+    subroutine s_compute_fd_divergence(div, fields, ix_s, iy_s, iz_s)
+
+        type(scalar_field), intent(INOUT) :: div
+        type(scalar_field), intent(IN) :: fields(1:3)
+        type(int_bounds_info), intent(IN) :: ix_s, iy_s, iz_s
+
+        integer :: x, y, z !< Generic loop iterators
+
+        real(kind(0d0)) :: divergence
+
+        !$acc parallel loop collapse(3) private(divergence)
+        do x = ix_s%beg, ix_s%end
+            do y = iy_s%beg, iy_s%end
+                do z = iz_s%beg, iz_s%end
+
+                    if (x == ix_s%beg) then
+                        divergence = (-3d0*fields(1)%sf(x, y, z) + 4d0*fields(1)%sf(x + 1, y, z) - fields(1)%sf(x + 2, y, z))/(x_cc(x + 2) - x_cc(x))
+                    else if (x == ix_s%end) then
+                        divergence = (+3d0*fields(1)%sf(x, y, z) - 4d0*fields(1)%sf(x - 1, y, z) + fields(1)%sf(x - 2, y, z))/(x_cc(x) - x_cc(x - 2))
+                    else
+                        divergence = (fields(1)%sf(x + 1, y, z) - fields(1)%sf(x - 1, y, z))/(x_cc(x + 1) - x_cc(x - 1))
+                    end if
+
+                    if (n > 0) then
+                        if (y == iy_s%beg) then
+                            divergence = divergence + (-3d0*fields(2)%sf(x, y, z) + 4d0*fields(2)%sf(x, y + 1, z) - fields(2)%sf(x, y + 2, z))/(y_cc(y + 2) - y_cc(y))
+                        else if (y == iy_s%end) then
+                            divergence = divergence + (+3d0*fields(2)%sf(x, y, z) - 4d0*fields(2)%sf(x, y - 1, z) + fields(2)%sf(x, y - 2, z))/(y_cc(y) - y_cc(y - 2))
+                        else
+                            divergence = divergence + (fields(2)%sf(x, y + 1, z) - fields(2)%sf(x, y - 1, z))/(y_cc(y + 1) - y_cc(y - 1))
+                        end if
+                    end if
+
+                    if (p > 0) then
+                        if (z == iz_s%beg) then
+                            divergence = divergence + (-3d0*fields(3)%sf(x, y, z) + 4d0*fields(3)%sf(x, y, z + 1) - fields(3)%sf(x, y, z + 2))/(z_cc(z + 2) - z_cc(z))
+                        else if (z == iz_s%end) then
+                            divergence = divergence + (+3d0*fields(3)%sf(x, y, z) - 4d0*fields(3)%sf(x, y, z - 1) + fields(2)%sf(x, y, z - 2))/(z_cc(z) - z_cc(z - 2))
+                        else
+                            divergence = divergence + (fields(3)%sf(x, y, z + 1) - fields(3)%sf(x, y, z - 1))/(z_cc(z + 1) - z_cc(z - 1))
+                        end if
+                    end if
+
+                    div%sf(x, y, z) = div%sf(x, y, z) + divergence
+
+                end do
+            end do
+        end do
+
+    end subroutine s_compute_fd_divergence
+
+    !>  The purpose of this subroutine is to compute the finite-
+    !!      difference coefficients for the centered schemes utilized
+    !!      in computations of first order spatial derivatives in the
+    !!      s-coordinate direction. The s-coordinate direction refers
+    !!      to the x-, y- or z-coordinate direction, depending on the
+    !!      subroutine's inputs. Note that coefficients of up to 4th
+    !!      order accuracy are available.
+    !!  @param q Number of cells in the s-coordinate direction
+    !!  @param s_cc Locations of the cell-centers in the s-coordinate direction
+    !!  @param fd_coeff_s Finite-diff. coefficients in the s-coordinate direction
+    subroutine s_compute_finite_difference_coefficients(q, s_cc, fd_coeff_s, buff_size, &
+                                                        fd_number_in, fd_order_in, offset_s)
+
+        integer :: lB, lE !< loop bounds
+        integer, intent(IN) :: q
+        integer, intent(IN) :: buff_size, fd_number_in, fd_order_in
+        type(int_bounds_info), optional, intent(IN) :: offset_s
+        real(kind(0d0)), allocatable, dimension(:, :), intent(INOUT) :: fd_coeff_s
+
+        real(kind(0d0)), &
+            dimension(-buff_size:q + buff_size), &
+            intent(IN) :: s_cc
+
+        integer :: i !< Generic loop iterator
+
+        if (present(offset_s)) then
+            lB = -offset_s%beg
+            lE = q + offset_s%end
+        else
+            lB = 0
+            lE = q
+        end if
+
+        if (allocated(fd_coeff_s)) deallocate (fd_coeff_s)
+        allocate (fd_coeff_s(-fd_number_in:fd_number_in, lb:lE))
+
+        ! Computing the 1st order finite-difference coefficients
+        if (fd_order_in == 1) then
+            do i = lB, lE
+                fd_coeff_s(-1, i) = 0d0
+                fd_coeff_s(0, i) = -1d0/(s_cc(i + 1) - s_cc(i))
+                fd_coeff_s(1, i) = -fd_coeff_s(0, i)
+            end do
+
+            ! Computing the 2nd order finite-difference coefficients
+        elseif (fd_order_in == 2) then
+            do i = lB, lE
+                fd_coeff_s(-1, i) = -1d0/(s_cc(i + 1) - s_cc(i - 1))
+                fd_coeff_s(0, i) = 0d0
+                fd_coeff_s(1, i) = -fd_coeff_s(-1, i)
+            end do
+
+            ! Computing the 4th order finite-difference coefficients
+        else
+            do i = lB, lE
+                fd_coeff_s(-2, i) = 1d0/(s_cc(i - 2) - 8d0*s_cc(i - 1) - s_cc(i + 2) + 8d0*s_cc(i + 1))
+                fd_coeff_s(-1, i) = -8d0*fd_coeff_s(-2, i)
+                fd_coeff_s(0, i) = 0d0
+                fd_coeff_s(1, i) = -fd_coeff_s(-1, i)
+                fd_coeff_s(2, i) = -fd_coeff_s(-2, i)
+            end do
+
+        end if
+
+    end subroutine s_compute_finite_difference_coefficients ! --------------
+
+end module m_finite_differences
+

src/common/m_helper.fpp

@@ -20,8 +20,7 @@ module m_helper
     implicit none
 
     private; 
-    public :: s_compute_finite_difference_coefficients, &
-              s_comp_n_from_prim, &
+    public :: s_comp_n_from_prim, &
               s_comp_n_from_cons, &
               s_initialize_nonpoly, &
               s_simpson, &
@@ -41,72 +40,6 @@ module m_helper
 
 contains
 
-    !>  The purpose of this subroutine is to compute the finite-
-        !!      difference coefficients for the centered schemes utilized
-        !!      in computations of first order spatial derivatives in the
-        !!      s-coordinate direction. The s-coordinate direction refers
-        !!      to the x-, y- or z-coordinate direction, depending on the
-        !!      subroutine's inputs. Note that coefficients of up to 4th
-        !!      order accuracy are available.
-        !!  @param q Number of cells in the s-coordinate direction
-        !!  @param s_cc Locations of the cell-centers in the s-coordinate direction
-        !!  @param fd_coeff_s Finite-diff. coefficients in the s-coordinate direction
-    subroutine s_compute_finite_difference_coefficients(q, s_cc, fd_coeff_s, buff_size, &
-                                                        fd_number_in, fd_order_in, offset_s)
-
-        integer, intent(in) :: q
-        real(kind(0d0)), allocatable, dimension(:, :), intent(inout) :: fd_coeff_s
-        integer, intent(in) :: buff_size, fd_number_in, fd_order_in
-        type(int_bounds_info), optional, intent(in) :: offset_s
-
-        real(kind(0d0)), &
-            dimension(-buff_size:q + buff_size), &
-            intent(IN) :: s_cc
-
-        integer :: lB, lE !< loop bounds
-        integer :: i !< Generic loop iterator
-
-        if (present(offset_s)) then
-            lB = -offset_s%beg
-            lE = q + offset_s%end
-        else
-            lB = 0
-            lE = q
-        end if
-
-        if (allocated(fd_coeff_s)) deallocate (fd_coeff_s)
-        allocate (fd_coeff_s(-fd_number_in:fd_number_in, lb:lE))
-
-        ! Computing the 1st order finite-difference coefficients
-        if (fd_order_in == 1) then
-            do i = lB, lE
-                fd_coeff_s(-1, i) = 0d0
-                fd_coeff_s(0, i) = -1d0/(s_cc(i + 1) - s_cc(i))
-                fd_coeff_s(1, i) = -fd_coeff_s(0, i)
-            end do
-
-            ! Computing the 2nd order finite-difference coefficients
-        elseif (fd_order_in == 2) then
-            do i = lB, lE
-                fd_coeff_s(-1, i) = -1d0/(s_cc(i + 1) - s_cc(i - 1))
-                fd_coeff_s(0, i) = 0d0
-                fd_coeff_s(1, i) = -fd_coeff_s(-1, i)
-            end do
-
-            ! Computing the 4th order finite-difference coefficients
-        else
-            do i = lB, lE
-                fd_coeff_s(-2, i) = 1d0/(s_cc(i - 2) - 8d0*s_cc(i - 1) - s_cc(i + 2) + 8d0*s_cc(i + 1))
-                fd_coeff_s(-1, i) = -8d0*fd_coeff_s(-2, i)
-                fd_coeff_s(0, i) = 0d0
-                fd_coeff_s(1, i) = -fd_coeff_s(-1, i)
-                fd_coeff_s(2, i) = -fd_coeff_s(-2, i)
-            end do
-
-        end if
-
-    end subroutine s_compute_finite_difference_coefficients
-
     !> Computes the bubble number density n from the primitive variables
         !! @param vftmp is the void fraction
         !! @param Rtmp is the  bubble radii

src/common/m_thermochem.fpp

@@ -0,0 +1,12 @@
+#:include 'case.fpp'
+
+module m_thermochem
+
+    #:if chemistry
+        use m_pyrometheus
+    #:else
+        integer, parameter :: num_species = 0
+        character(len=:), allocatable, dimension(:) :: species_names
+    #:endif
+
+end module m_thermochem