From 3b53b1bbff3a4c7563e7436d67f0ffd6cbae4ed2 Mon Sep 17 00:00:00 2001 From: Henry LE BERRE Date: Mon, 29 Jul 2024 18:44:29 +0200 Subject: [PATCH] Upstream CheMFC infrastructure --- .gitattributes | 1 + .github/workflows/docs.yml | 9 +- .github/workflows/test.yml | 4 +- CMakeLists.txt | 1 + README.md | 8 +- docs/documentation/getting-started.md | 8 +- misc/fpp_to_fypp.sh | 6 + src/common/include/macros.fpp | 18 +- src/common/m_derived_types.fpp | 18 +- src/common/m_finite_differences.fpp | 127 ++++++++++++ src/common/m_helper.fpp | 69 +------ src/common/m_thermochem.fpp | 12 ++ src/common/m_variables_conversion.fpp | 185 ++++++++++++----- src/post_process/m_global_parameters.fpp | 29 +++ src/post_process/m_start_up.f90 | 23 +++ ...n_variables.f90 => m_assign_variables.fpp} | 75 ++++++- src/pre_process/m_check_patches.fpp | 10 +- src/pre_process/m_checker.fpp | 9 + src/pre_process/m_data_output.fpp | 74 +++++-- src/pre_process/m_global_parameters.fpp | 28 +++ src/pre_process/m_patches.fpp | 15 +- src/simulation/m_acoustic_src.fpp | 2 + src/simulation/m_boundary_conditions.fpp | 7 +- src/simulation/m_chemistry.fpp | 188 ++++++++++++++++++ src/simulation/m_compute_cbc.fpp | 2 +- src/simulation/m_data_output.fpp | 54 ++++- src/simulation/m_derived_variables.f90 | 3 + src/simulation/m_global_parameters.fpp | 34 +++- src/simulation/m_mpi_proxy.fpp | 10 +- src/simulation/m_rhs.fpp | 55 ++++- src/simulation/m_riemann_solvers.fpp | 19 +- src/simulation/m_start_up.fpp | 23 ++- src/simulation/m_time_steppers.fpp | 20 ++ src/simulation/m_viscous.fpp | 7 +- toolchain/mfc/build.py | 3 + toolchain/mfc/case.py | 9 +- toolchain/mfc/run/case_dicts.py | 25 +-- toolchain/mfc/run/input.py | 34 ++++ toolchain/requirements.txt | 26 ++- 39 files changed, 1039 insertions(+), 211 deletions(-) create mode 100644 misc/fpp_to_fypp.sh create mode 100644 src/common/m_finite_differences.fpp create mode 100644 src/common/m_thermochem.fpp rename src/pre_process/{m_assign_variables.f90 => m_assign_variables.fpp} (90%) create mode 100644 src/simulation/m_chemistry.fpp diff --git a/.gitattributes b/.gitattributes index e0de9cc53..442e4b432 100644 --- a/.gitattributes +++ b/.gitattributes @@ -1 +1,2 @@ /tests/**/* linguist-generated=true +/toolchain/mechanisms/* linguist-generated=true diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index 6d5ae2180..cd329acf4 100644 --- a/.github/workflows/docs.yml +++ b/.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 diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 248c07211..de063affb 100644 --- a/.github/workflows/test.yml +++ b/.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 - diff --git a/CMakeLists.txt b/CMakeLists.txt index 7d05c9c19..8e5e8138e 100644 --- a/CMakeLists.txt +++ b/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}" diff --git a/README.md b/README.md index f816174f7..e337ae36d 100644 --- a/README.md +++ b/README.md @@ -58,7 +58,7 @@ 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 @@ -66,6 +66,12 @@ Put it to a convenient directory via 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) diff --git a/docs/documentation/getting-started.md b/docs/documentation/getting-started.md index a68ebe927..a47074e5f 100644 --- a/docs/documentation/getting-started.md +++ b/docs/documentation/getting-started.md @@ -101,10 +101,14 @@ You will also have access to the `.sln` Microsoft Visual Studio solution files f

MacOS

 -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. diff --git a/misc/fpp_to_fypp.sh b/misc/fpp_to_fypp.sh new file mode 100644 index 000000000..50bb797e3 --- /dev/null +++ b/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 diff --git a/src/common/include/macros.fpp b/src/common/include/macros.fpp index e44290c26..27a1bf338 100644 --- a/src/common/include/macros.fpp +++ b/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 diff --git a/src/common/m_derived_types.fpp b/src/common/m_derived_types.fpp index d8c067450..65d409474 100644 --- a/src/common/m_derived_types.fpp +++ b/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 diff --git a/src/common/m_finite_differences.fpp b/src/common/m_finite_differences.fpp new file mode 100644 index 000000000..9eb65a121 --- /dev/null +++ b/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 + diff --git a/src/common/m_helper.fpp b/src/common/m_helper.fpp index f49b4adc7..fc26005fd 100644 --- a/src/common/m_helper.fpp +++ b/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 diff --git a/src/common/m_thermochem.fpp b/src/common/m_thermochem.fpp new file mode 100644 index 000000000..87baf6f8d --- /dev/null +++ b/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 diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index 180abefd6..35c1ad975 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -21,6 +21,9 @@ module m_variables_conversion use m_helper_basic !< Functions to compare floating point numbers use m_helper + + use m_thermochem + ! ========================================================================== implicit none @@ -124,7 +127,7 @@ contains !! @param pres Pressure to calculate !! @param stress Shear Stress !! @param mom Momentum - subroutine s_compute_pressure(energy, alf, dyn_p, pi_inf, gamma, rho, qv, pres, stress, mom, G) + subroutine s_compute_pressure(energy, alf, dyn_p, pi_inf, gamma, rho, qv, rhoYks, pres, stress, mom, G) !$acc routine seq real(kind(0d0)), intent(in) :: energy, alf @@ -133,46 +136,67 @@ contains real(kind(0d0)), intent(out) :: pres real(kind(0d0)), intent(in), optional :: stress, mom, G + ! Chemistry + integer :: i + real(kind(0d0)), dimension(1:num_species), intent(in) :: rhoYks real(kind(0d0)) :: E_e + real(kind(0d0)) :: T + real(kind(0d0)), dimension(1:num_species) :: Y_rs integer :: s !< Generic loop iterator - ! Depending on model_eqns and bubbles, the appropriate procedure - ! for computing pressure is targeted by the procedure pointer + #:if not chemistry + ! 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 - pres = (energy - dyn_p - pi_inf - qv)/gamma - else if ((model_eqns /= 4) .and. bubbles) then - pres = ((energy - dyn_p)/(1.d0 - alf) - pi_inf - qv)/gamma - else - pres = (pref + pi_inf)* & - (energy/ & - (rhoref*(1 - alf)) & - )**(1/gamma + 1) - pi_inf - end if + if ((model_eqns /= 4) .and. (bubbles .neqv. .true.)) then + pres = (energy - dyn_p - pi_inf - qv)/gamma + else if ((model_eqns /= 4) .and. bubbles) then + pres = ((energy - dyn_p)/(1.d0 - alf) - pi_inf - qv)/gamma + else + pres = (pref + pi_inf)* & + (energy/ & + (rhoref*(1 - alf)) & + )**(1/gamma + 1) - pi_inf + end if - if (hypoelasticity .and. present(G)) then - ! calculate elastic contribution to Energy - E_e = 0d0 - do s = stress_idx%beg, stress_idx%end - if (G > 0) then - E_e = E_e + ((stress/rho)**2d0)/(4d0*G) - ! Additional terms in 2D and 3D - if ((s == stress_idx%beg + 1) .or. & - (s == stress_idx%beg + 3) .or. & - (s == stress_idx%beg + 4)) then + if (hypoelasticity .and. present(G)) then + ! calculate elastic contribution to Energy + E_e = 0d0 + do s = stress_idx%beg, stress_idx%end + if (G > 0) then E_e = E_e + ((stress/rho)**2d0)/(4d0*G) + ! Additional terms in 2D and 3D + if ((s == stress_idx%beg + 1) .or. & + (s == stress_idx%beg + 3) .or. & + (s == stress_idx%beg + 4)) then + E_e = E_e + ((stress/rho)**2d0)/(4d0*G) + end if end if - end if + end do + + pres = ( & + energy - & + 0.5d0*(mom**2.d0)/rho - & + pi_inf - qv - E_e & + )/gamma + + end if + + #:else + !$acc loop seq + do i = 1, num_species + Y_rs(i) = rhoYks(i)/rho end do - pres = ( & - energy - & - 0.5d0*(mom**2.d0)/rho - & - pi_inf - qv - E_e & - )/gamma + if (sum(Y_rs) > 1d-16) then + call get_temperature(.true., energy - dyn_p, 1200d0, Y_rs, T) + call get_pressure(rho, T, Y_rs, pres) + else + pres = 0d0 + end if - end if + #:endif end subroutine s_compute_pressure @@ -856,11 +880,13 @@ contains real(kind(0d0)), dimension(:), allocatable :: nRtmp #:endif + real(kind(0d0)) :: rhoYks(1:num_species) + real(kind(0d0)) :: vftmp, nR3, nbub_sc, R3tmp real(kind(0d0)) :: G_K - real(kind(0d0)) :: pres + real(kind(0d0)) :: pres, Yksum integer :: i, j, k, l, q !< Generic loop iterators @@ -882,7 +908,7 @@ contains end if #:endif - !$acc parallel loop collapse(3) gang vector default(present) private(alpha_K, alpha_rho_K, Re_K, nRtmp, rho_K, gamma_K, pi_inf_K, qv_K, dyn_pres_K, R3tmp) + !$acc parallel loop collapse(3) gang vector default(present) private(alpha_K, alpha_rho_K, Re_K, nRtmp, rho_K, gamma_K, pi_inf_K, qv_K, dyn_pres_K, R3tmp, rhoyks) do l = izb, ize do k = iyb, iye do j = ixb, ixe @@ -894,11 +920,6 @@ contains alpha_K(i) = qK_cons_vf(advxb + i - 1)%sf(j, k, l) end do - !$acc loop seq - do i = 1, contxe - qK_prim_vf(i)%sf(j, k, l) = qK_cons_vf(i)%sf(j, k, l) - end do - if (model_eqns /= 4) then #ifdef MFC_SIMULATION ! If in simulation, use acc mixture subroutines @@ -924,6 +945,39 @@ contains #endif end if + if (chemistry) then + rho_K = 0d0 + !$acc loop seq + do i = chemxb, chemxe + !print*, j,k,l, qK_cons_vf(i)%sf(j, k, l) + rho_K = rho_K + max(0d0, qK_cons_vf(i)%sf(j, k, l)) + end do + + !$acc loop seq + do i = 1, contxe + qK_prim_vf(i)%sf(j, k, l) = rho_K + end do + + Yksum = 0d0 + !$acc loop seq + do i = chemxb, chemxe + qK_prim_vf(i)%sf(j, k, l) = max(0d0, qK_cons_vf(i)%sf(j, k, l)/rho_K) + Yksum = Yksum + qK_prim_vf(i)%sf(j, k, l) + end do + + !$acc loop seq + do i = chemxb, chemxe + qK_prim_vf(i)%sf(j, k, l) = qK_prim_vf(i)%sf(j, k, l)/Yksum + end do + + qK_prim_vf(tempxb)%sf(j, k, l) = qK_cons_vf(tempxb)%sf(j, k, l) + else + !$acc loop seq + do i = 1, contxe + qK_prim_vf(i)%sf(j, k, l) = qK_cons_vf(i)%sf(j, k, l) + end do + end if + #ifdef MFC_SIMULATION rho_K = max(rho_K, sgm_eps) #endif @@ -941,9 +995,16 @@ contains end if end do + if (chemistry) then + !$acc loop seq + do i = 1, num_species + rhoYks(i) = qK_cons_vf(chemxb + i - 1)%sf(j, k, l) + end do + end if + call s_compute_pressure(qK_cons_vf(E_idx)%sf(j, k, l), & qK_cons_vf(alf_idx)%sf(j, k, l), & - dyn_pres_K, pi_inf_K, gamma_K, rho_K, qv_K, pres) + dyn_pres_K, pi_inf_K, gamma_K, rho_K, qv_K, rhoYks, pres) qK_prim_vf(E_idx)%sf(j, k, l) = pres @@ -1053,6 +1114,10 @@ contains real(kind(0d0)), dimension(2) :: Re_K integer :: i, j, k, l, q !< Generic loop iterators + integer :: spec + + real(kind(0d0)), dimension(num_species) :: Ys + real(kind(0d0)) :: temperature, e_mix, mix_mol_weight, T #ifndef MFC_SIMULATION ! Converting the primitive variables to the conservative variables @@ -1086,21 +1151,37 @@ contains q_prim_vf(i)%sf(j, k, l)/2d0 end do - ! Computing the energy from the pressure - if ((model_eqns /= 4) .and. (bubbles .neqv. .true.)) then - ! E = Gamma*P + \rho u u /2 + \pi_inf + (\alpha\rho qv) + #:if chemistry + do i = chemxb, chemxe + Ys(i - chemxb + 1) = q_prim_vf(i)%sf(j, k, l) + q_cons_vf(i)%sf(j, k, l) = rho*q_prim_vf(i)%sf(j, k, l) + end do + + call get_mixture_molecular_weight(Ys, mix_mol_weight) + T = q_prim_vf(E_idx)%sf(j, k, l)*mix_mol_weight/(gas_constant*rho) + call get_mixture_energy_mass(T, Ys, e_mix) + q_cons_vf(E_idx)%sf(j, k, l) = & - gamma*q_prim_vf(E_idx)%sf(j, k, l) + dyn_pres + pi_inf & - + qv - else if ((model_eqns /= 4) .and. (bubbles)) then - ! \tilde{E} = dyn_pres + (1-\alf)(\Gamma p_l + \Pi_inf) - q_cons_vf(E_idx)%sf(j, k, l) = dyn_pres + & - (1.d0 - q_prim_vf(alf_idx)%sf(j, k, l))* & - (gamma*q_prim_vf(E_idx)%sf(j, k, l) + pi_inf) - else - !Tait EOS, no conserved energy variable - q_cons_vf(E_idx)%sf(j, k, l) = 0. - end if + dyn_pres + e_mix + + q_cons_vf(tempxb)%sf(j, k, l) = T + #:else + ! Computing the energy from the pressure + if ((model_eqns /= 4) .and. (bubbles .neqv. .true.)) then + ! E = Gamma*P + \rho u u /2 + \pi_inf + (\alpha\rho qv) + q_cons_vf(E_idx)%sf(j, k, l) = & + gamma*q_prim_vf(E_idx)%sf(j, k, l) + dyn_pres + pi_inf & + + qv + else if ((model_eqns /= 4) .and. (bubbles)) then + ! \tilde{E} = dyn_pres + (1-\alf)(\Gamma p_l + \Pi_inf) + q_cons_vf(E_idx)%sf(j, k, l) = dyn_pres + & + (1.d0 - q_prim_vf(alf_idx)%sf(j, k, l))* & + (gamma*q_prim_vf(E_idx)%sf(j, k, l) + pi_inf) + else + !Tait EOS, no conserved energy variable + q_cons_vf(E_idx)%sf(j, k, l) = 0. + end if + #:endif ! Computing the internal energies from the pressure and continuities if (model_eqns == 3) then diff --git a/src/post_process/m_global_parameters.fpp b/src/post_process/m_global_parameters.fpp index 8a6b2805f..312759ab2 100644 --- a/src/post_process/m_global_parameters.fpp +++ b/src/post_process/m_global_parameters.fpp @@ -2,6 +2,8 @@ !! @file m_global_parameters.f90 !! @brief Contains module m_global_parameters +#:include 'case.fpp' + !> @brief This module contains all of the parameters characterizing the !! computational domain, simulation algorithm, stiffened equation of !! state and finally, the formatted database file(s) structure. @@ -15,6 +17,9 @@ module m_global_parameters use m_derived_types !< Definitions of the derived types use m_helper_basic !< Functions to compare floating point numbers + + use m_thermochem !< Thermodynamic and chemical properties module + ! ========================================================================== implicit none @@ -104,6 +109,7 @@ module m_global_parameters logical :: mixture_err !< Mixture error limiter logical :: alt_soundspeed !< Alternate sound speed logical :: hypoelasticity !< Turn hypoelasticity on + logical, parameter :: chemistry = .${chemistry}$. !< Chemistry modeling !> @} !> @name Annotations of the structure, i.e. the organization, of the state vectors @@ -120,6 +126,8 @@ module m_global_parameters integer :: pi_inf_idx !< Index of liquid stiffness func. eqn. type(int_bounds_info) :: stress_idx !< Indices of elastic stresses integer :: c_idx !< Index of color function + type(int_bounds_info) :: chemistry_idx !< Indexes of first & last concentration eqns. + type(int_bounds_info) :: temperature_idx !< Indexes of first & last temperature eqns. !> @} !> @name Boundary conditions in the x-, y- and z-coordinate directions @@ -204,6 +212,8 @@ module m_global_parameters logical :: schlieren_wrt logical :: cf_wrt logical :: ib + logical :: chem_wrt_Y(1:num_species) + logical :: chem_wrt_T !> @} real(kind(0d0)), dimension(num_fluids_max) :: schlieren_alpha !< @@ -265,6 +275,8 @@ module m_global_parameters integer :: intxb, intxe integer :: bubxb, bubxe integer :: strxb, strxe + integer :: chemxb, chemxe + integer :: tempxb, tempxe !> @} contains @@ -336,6 +348,8 @@ contains rho_wrt = .false. mom_wrt = .false. vel_wrt = .false. + chem_wrt_Y = .false. + chem_wrt_T = .false. flux_lim = dflt_int flux_wrt = .false. parallel_io = .false. @@ -596,6 +610,16 @@ contains end if end if + if (chemistry) then + chemistry_idx%beg = sys_size + 1 + chemistry_idx%end = sys_size + num_species + sys_size = chemistry_idx%end + + temperature_idx%beg = sys_size + 1 + temperature_idx%end = sys_size + 1 + sys_size = temperature_idx%end + end if + momxb = mom_idx%beg momxe = mom_idx%end advxb = adv_idx%beg @@ -608,6 +632,11 @@ contains strxe = stress_idx%end intxb = internalEnergies_idx%beg intxe = internalEnergies_idx%end + chemxb = chemistry_idx%beg + chemxe = chemistry_idx%end + tempxb = temperature_idx%beg + tempxe = temperature_idx%end + ! ================================================================== #ifdef MFC_MPI diff --git a/src/post_process/m_start_up.f90 b/src/post_process/m_start_up.f90 index 2dd94da62..78b42a7bf 100644 --- a/src/post_process/m_start_up.f90 +++ b/src/post_process/m_start_up.f90 @@ -34,6 +34,12 @@ module m_start_up use m_checker_common use m_checker + + use m_thermochem !< Procedures used to compute thermodynamic + !! quantities + + use m_finite_differences + ! ========================================================================== implicit none @@ -64,6 +70,7 @@ subroutine s_read_input_file weno_order, bc_x, & bc_y, bc_z, fluid_pp, format, precision, & hypoelasticity, G, & + chem_wrt_Y, chem_wrt_T, & alpha_rho_wrt, rho_wrt, mom_wrt, vel_wrt, & E_wrt, pres_wrt, alpha_wrt, gamma_wrt, & heat_ratio_wrt, pi_inf_wrt, pres_inf_wrt, & @@ -174,6 +181,7 @@ subroutine s_perform_time_step(t_step) ! Converting the conservative variables to the primitive ones call s_convert_conservative_to_primitive_variables(q_cons_vf, q_prim_vf) + end subroutine s_perform_time_step subroutine s_save_data(t_step, varname, pres, c, H) @@ -287,6 +295,21 @@ subroutine s_save_data(t_step, varname, pres, c, H) end do ! ---------------------------------------------------------------------- + ! Adding the species' concentrations to the formatted database file ---- + do i = 1, num_species + if (chem_wrt_Y(i) .or. prim_vars_wrt) then + q_sf = q_prim_vf(chemxb + i - 1)%sf(-offset_x%beg:m + offset_x%end, & + -offset_y%beg:n + offset_y%end, & + -offset_z%beg:p + offset_z%end) + + write (varname, '(A,A)') 'Y_', trim(species_names(i)) + call s_write_variable_to_formatted_database_file(varname, t_step) + + varname(:) = ' ' + + end if + end do + ! Adding the flux limiter function to the formatted database file do i = 1, E_idx - mom_idx%beg if (flux_wrt(i)) then diff --git a/src/pre_process/m_assign_variables.f90 b/src/pre_process/m_assign_variables.fpp similarity index 90% rename from src/pre_process/m_assign_variables.f90 rename to src/pre_process/m_assign_variables.fpp index c544ffa2e..5aa19e6a3 100644 --- a/src/pre_process/m_assign_variables.f90 +++ b/src/pre_process/m_assign_variables.fpp @@ -1,6 +1,9 @@ !> !! @file m_assign_variables.f90 !! @brief Contains module m_assign_variables + +#:include 'case.fpp' + module m_assign_variables ! Dependencies ============================================================= @@ -10,7 +13,9 @@ module m_assign_variables use m_variables_conversion ! Subroutines to change the state variables from - use m_helper_basic !< Functions to compare floating point numbers + use m_helper_basic !< Functions to compare floating point numbers + + use m_thermochem !< Thermodynamic and chemical properties ! one form to another ! ========================================================================== @@ -116,6 +121,8 @@ subroutine s_assign_patch_mixture_primitive_variables(patch_id, j, k, l, & real(kind(0d0)) :: gamma !< specific heat ratio function real(kind(0d0)) :: x_centroid, y_centroid real(kind(0d0)) :: epsilon, beta + real(kind(0d0)) :: Ys(1:num_species) + real(kind(0d0)) :: mean_molecular_weight integer :: smooth_patch_id integer :: i !< generic loop operator @@ -158,6 +165,37 @@ subroutine s_assign_patch_mixture_primitive_variables(patch_id, j, k, l, & eta*patch_icpp(patch_id)%pi_inf & + (1d0 - eta)*patch_icpp(smooth_patch_id)%pi_inf + ! Species Concentrations + #:if chemistry + block + real(kind(0d0)) :: sum, term + + ! Accumulating the species concentrations + sum = 0d0 + do i = 1, num_species + term = & + eta*patch_icpp(patch_id)%Y(i) & + + (1d0 - eta)*patch_icpp(smooth_patch_id)%Y(i) + q_prim_vf(chemxb + i - 1)%sf(j, k, l) = term + sum = sum + term + end do + + sum = max(sum, verysmall) + + ! Normalizing the species concentrations + do i = 1, num_species + q_prim_vf(chemxb + i - 1)%sf(j, k, l) = & + q_prim_vf(chemxb + i - 1)%sf(j, k, l)/sum + Ys(i) = q_prim_vf(chemxb + i - 1)%sf(j, k, l) + end do + end block + + call get_mixture_molecular_weight(Ys, mean_molecular_weight) + q_prim_vf(tempxb)%sf(j, k, l) = & + q_prim_vf(E_idx)%sf(j, k, l)*mean_molecular_weight & + /(gas_constant*q_prim_vf(1)%sf(j, k, l)) + #:endif + ! Updating the patch identities bookkeeping variable if (1d0 - eta < 1d-16) patch_id_fp(j, k, l) = patch_id @@ -282,6 +320,9 @@ subroutine s_assign_patch_species_primitive_variables(patch_id, j, k, l, & real(kind(0d0)) :: x_centroid, y_centroid real(kind(0d0)) :: epsilon, beta + real(kind(0d0)) :: Ys(1:num_species) + real(kind(0d0)) :: mean_molecular_weight + real(kind(0d0)), dimension(sys_size) :: orig_prim_vf !< !! Vector to hold original values of cell for smoothing purposes @@ -501,6 +542,38 @@ subroutine s_assign_patch_species_primitive_variables(patch_id, j, k, l, & + (1d0 - eta)*orig_prim_vf(i + cont_idx%end)) end do + ! Species Concentrations + #:if chemistry + block + real(kind(0d0)) :: sum, term + + ! Accumulating the species concentrations + sum = 0d0 + do i = 1, num_species + term = & + eta*patch_icpp(patch_id)%Y(i) & + + (1d0 - eta)*patch_icpp(smooth_patch_id)%Y(i) + q_prim_vf(chemxb + i - 1)%sf(j, k, l) = term + sum = sum + term + end do + + if (sum < verysmall) then + sum = 1d0 + end if + + ! Normalizing the species concentrations + do i = 1, num_species + q_prim_vf(chemxb + i - 1)%sf(j, k, l) = & + q_prim_vf(chemxb + i - 1)%sf(j, k, l)/sum + Ys(i) = q_prim_vf(chemxb + i - 1)%sf(j, k, l) + end do + end block + + call get_mixture_molecular_weight(Ys, mean_molecular_weight) + q_prim_vf(tempxb)%sf(j, k, l) = & + q_prim_vf(E_idx)%sf(j, k, l)*mean_molecular_weight/(gas_constant*q_prim_vf(1)%sf(j, k, l)) + #:endif + ! Set streamwise velocity to hyperbolic tangent function of y if (mixlayer_vel_profile) then q_prim_vf(1 + cont_idx%end)%sf(j, k, l) = & diff --git a/src/pre_process/m_check_patches.fpp b/src/pre_process/m_check_patches.fpp index 6e1e69396..a9ddc3883 100644 --- a/src/pre_process/m_check_patches.fpp +++ b/src/pre_process/m_check_patches.fpp @@ -1,3 +1,5 @@ +#:include 'macros.fpp' + !> @brief This module contains subroutines that read, and check consistency !! of, the user provided inputs and data. @@ -505,6 +507,7 @@ contains subroutine s_check_active_patch_primitive_variables(patch_id) integer, intent(in) :: patch_id + call s_int_to_str(patch_id, iStr) @:PROHIBIT(f_is_default(patch_icpp(patch_id)%vel(1)), & @@ -537,7 +540,12 @@ contains "Patch "//trim(iStr)//": alpha(num_fluids) must be set") end if - end subroutine s_check_active_patch_primitive_variables + if (chemistry) then + !@:ASSERT(all(patch_icpp(patch_id)%Y(1:num_species) >= 0d0), "Patch " // trim(iStr) // ".") + !@:ASSERT(any(patch_icpp(patch_id)%Y(1:num_species) > verysmall), "Patch " // trim(iStr) // ".") + end if + + end subroutine s_check_active_patch_primitive_variables ! -------------- !> This subroutine verifies that the primitive variables !! associated with the given inactive patch remain unaltered diff --git a/src/pre_process/m_checker.fpp b/src/pre_process/m_checker.fpp index fbc15dc8b..efd1ef190 100644 --- a/src/pre_process/m_checker.fpp +++ b/src/pre_process/m_checker.fpp @@ -30,6 +30,7 @@ contains call s_check_inputs_grid_stretching call s_check_inputs_qbmm_and_polydisperse call s_check_inputs_perturb_density + call s_check_inputs_chemistry call s_check_inputs_misc end subroutine s_check_inputs @@ -170,6 +171,14 @@ contains end do end subroutine s_check_inputs_perturb_density + subroutine s_check_inputs_chemistry + + if (chemistry) then + @:ASSERT(num_species > 0) + end if + + end subroutine s_check_inputs_chemistry + !> Checks miscellaneous constraints !! (mixlayer_vel_profile and mixlayer_perturb) subroutine s_check_inputs_misc diff --git a/src/pre_process/m_data_output.fpp b/src/pre_process/m_data_output.fpp index 41261eb34..e17773787 100644 --- a/src/pre_process/m_data_output.fpp +++ b/src/pre_process/m_data_output.fpp @@ -97,7 +97,7 @@ contains character(LEN=len_trim(t_step_dir) + name_len) :: file_loc !< !! Generic string used to store the address of a particular file - integer :: i, j, k, l, r !< Generic loop iterator + integer :: i, j, k, l, r, c !< Generic loop iterator integer :: t_step real(kind(0d0)), dimension(nb) :: nRtmp !< Temporary bubble concentration @@ -109,6 +109,8 @@ contains real(kind(0d0)) :: nR3 real(kind(0d0)) :: ntmp + real(kind(0d0)) :: rhoYks(1:num_species) !< Temporary species mass fractions + t_step = 0 ! Outputting the Locations of the Cell-boundaries ================== @@ -240,14 +242,27 @@ contains open (2, FILE=trim(file_loc)) do j = 0, m + + if (chemistry) then + do c = 1, num_species + rhoYks(c) = q_cons_vf(chemxb + c - 1)%sf(j, 0, 0) + end do + end if + call s_convert_to_mixture_variables(q_cons_vf, j, 0, 0, rho, gamma, pi_inf, qv) lit_gamma = 1d0/gamma + 1d0 - if (((i >= cont_idx%beg) .and. (i <= cont_idx%end)) & - .or. & - ((i >= adv_idx%beg) .and. (i <= adv_idx%end)) & - ) then + if ((i >= chemxb) .and. (i <= chemxe)) then + write (2, FMT) x_cb(j), q_cons_vf(i)%sf(j, 0, 0)/rho + else if (((i >= cont_idx%beg) .and. (i <= cont_idx%end)) & + .or. & + ((i >= adv_idx%beg) .and. (i <= adv_idx%end)) & + .or. & + ((i >= chemxb) .and. (i <= chemxe)) & + .or. & + ((i == tempxb)) & + ) then write (2, FMT) x_cb(j), q_cons_vf(i)%sf(j, 0, 0) else if (i == mom_idx%beg) then !u write (2, FMT) x_cb(j), q_cons_vf(mom_idx%beg)%sf(j, 0, 0)/rho @@ -258,7 +273,7 @@ contains q_cons_vf(E_idx)%sf(j, 0, 0), & q_cons_vf(alf_idx)%sf(j, 0, 0), & 0.5d0*(q_cons_vf(mom_idx%beg)%sf(j, 0, 0)**2.d0)/rho, & - pi_inf, gamma, rho, qv, pres) + pi_inf, gamma, rho, qv, rhoYks, pres) write (2, FMT) x_cb(j), pres else if ((i >= bub_idx%beg) .and. (i <= bub_idx%end) .and. bubbles) then @@ -735,9 +750,12 @@ contains subroutine s_initialize_data_output_module ! Generic string used to store the address of a particular file character(LEN=len_trim(case_dir) + 2*name_len) :: file_loc + character(len=15) :: temp + character(LEN=1), dimension(3) :: coord = (/'x', 'y', 'z'/) ! Generic logical used to check the existence of directories logical :: dir_check + integer :: i if (parallel_io .neqv. .true.) then ! Setting the address of the time-step directory @@ -781,15 +799,43 @@ contains end if - open (1, FILE='equations.dat', STATUS='unknown') + open (1, FILE='indices.dat', STATUS='unknown') + + write (1, '(A)') "Warning: The creation of file is currently experimental." + write (1, '(A)') "This file may contain errors and not support all features." + + write (1, '(A3,A20,A20)') "#", "Conservative", "Primitive" + write (1, '(A)') "-------------------------------------------" + do i = contxb, contxe + write (temp, '(I0)') i - contxb + 1 + write (1, '(I3,A20,A20)') i, "\alpha_{"//trim(temp)//"} \rho_{"//trim(temp)//"}", "\alpha_{"//trim(temp)//"} \rho" + end do + do i = momxb, momxe + write (1, '(I3,A20,A20)') i, "\rho u_"//coord(i - momxb + 1), "u_"//coord(i - momxb + 1) + end do + do i = E_idx, E_idx + write (1, '(I3,A20,A20)') i, "\rho U", "p" + end do + do i = advxb, advxe + write (temp, '(I0)') i - contxb + 1 + write (1, '(I3,A20,A20)') i, "\alpha_{"//trim(temp)//"}", "\alpha_{"//trim(temp)//"}" + end do + if (chemistry) then + do i = 1, num_species + write (1, '(I3,A20,A20)') chemxb + i - 1, "Y_{"//trim(species_names(i))//"} \rho", "Y_{"//trim(species_names(i))//"}" + end do + end if - write (1, '(A)') "Equations: " - if (momxb /= 0) write (1, '(A,I3,I3)') " * Momentum: ", momxb, momxe - if (advxb /= 0) write (1, '(A,I3,I3)') " * Advection: ", advxb, advxe - if (contxb /= 0) write (1, '(A,I3,I3)') " * Continuity: ", contxb, contxe - if (bubxb /= 0) write (1, '(A,I3,I3)') " * Bubbles: ", bubxb, bubxe - if (strxb /= 0) write (1, '(A,I3,I3)') " * Stress: ", strxb, strxe - if (intxb /= 0) write (1, '(A,I3,I3)') " * Internal Energies: ", intxb, intxe + write (1, '(A)') "" + if (momxb /= 0) write (1, '("[",I2,",",I2,"]",A)') momxb, momxe, " Momentum" + if (E_idx /= 0) write (1, '("[",I2,",",I2,"]",A)') E_idx, E_idx, " Energy/Pressure" + if (advxb /= 0) write (1, '("[",I2,",",I2,"]",A)') advxb, advxe, " Advection" + if (contxb /= 0) write (1, '("[",I2,",",I2,"]",A)') contxb, contxe, " Continuity" + if (bubxb /= 0) write (1, '("[",I2,",",I2,"]",A)') bubxb, bubxe, " Bubbles" + if (strxb /= 0) write (1, '("[",I2,",",I2,"]",A)') strxb, strxe, " Stress" + if (intxb /= 0) write (1, '("[",I2,",",I2,"]",A)') intxb, intxe, " Internal Energies" + if (chemxb /= 0) write (1, '("[",I2,",",I2,"]",A)') chemxb, chemxe, " Chemistry" + if (tempxb /= 0) write (1, '("[",I2,",",I2,"]",A)') tempxb, tempxe, " Temperature" close (1) diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index e517a5f46..6f15347bd 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/src/pre_process/m_global_parameters.fpp @@ -2,6 +2,8 @@ !! @file m_global_parameters.f90 !! @brief Contains module m_global_parameters +#:include 'case.fpp' + !> @brief This module contains all of the parameters characterizing the !! computational domain, simulation algorithm, initial condition !! and the stiffened equation of state. @@ -15,6 +17,9 @@ module m_global_parameters use m_derived_types ! Definitions of the derived types use m_helper_basic ! Functions to compare floating point numbers + + use m_thermochem ! Thermodynamic and chemical properties + ! ========================================================================== implicit none @@ -85,6 +90,7 @@ module m_global_parameters integer :: sys_size !< Number of unknowns in the system of equations integer :: weno_order !< Order of accuracy for the WENO reconstruction logical :: hypoelasticity !< activate hypoelasticity + logical, parameter :: chemistry = .${chemistry}$. !< Chemistry modeling ! Annotations of the structure, i.e. the organization, of the state vectors type(int_bounds_info) :: cont_idx !< Indexes of first & last continuity eqns. @@ -99,6 +105,8 @@ module m_global_parameters integer :: pi_inf_idx !< Index of liquid stiffness func. eqn. type(int_bounds_info) :: stress_idx !< Indexes of elastic shear stress eqns. integer :: c_idx !< Index of the color function + type(int_bounds_info) :: chemistry_idx !< Indexes of first & last concentration eqns. + type(int_bounds_info) :: temperature_idx !< Indexes of first & last temperature eqns. type(int_bounds_info) :: bc_x, bc_y, bc_z !< !! Boundary conditions in the x-, y- and z-coordinate directions @@ -221,6 +229,8 @@ module m_global_parameters integer :: intxb, intxe integer :: bubxb, bubxe integer :: strxb, strxe + integer :: chemxb, chemxe + integer :: tempxb, tempxe !> @} integer, allocatable, dimension(:, :, :) :: logic_grid @@ -359,6 +369,10 @@ contains patch_icpp(i)%m0 = dflt_real patch_icpp(i)%hcid = dflt_int + + if (chemistry) then + patch_icpp(i)%Y(:) = 0d0 + end if end do ! Tait EOS @@ -681,6 +695,16 @@ contains end if end if + if (chemistry) then + chemistry_idx%beg = sys_size + 1 + chemistry_idx%end = sys_size + num_species + sys_size = chemistry_idx%end + + temperature_idx%beg = sys_size + 1 + temperature_idx%end = sys_size + 1 + sys_size = temperature_idx%end + end if + momxb = mom_idx%beg momxe = mom_idx%end advxb = adv_idx%beg @@ -693,6 +717,10 @@ contains strxe = stress_idx%end intxb = internalEnergies_idx%beg intxe = internalEnergies_idx%end + chemxb = chemistry_idx%beg + chemxe = chemistry_idx%end + tempxb = temperature_idx%beg + tempxe = temperature_idx%end ! ================================================================== diff --git a/src/pre_process/m_patches.fpp b/src/pre_process/m_patches.fpp index 0fe993e8a..5baf28cc8 100644 --- a/src/pre_process/m_patches.fpp +++ b/src/pre_process/m_patches.fpp @@ -988,7 +988,20 @@ contains y_boundary%end >= y_cc(j)) & then - patch_id_fp(i, j, 0) = patch_id + call s_assign_patch_primitive_variables(patch_id, i, j, 0, & + eta, q_prim_vf, patch_id_fp) + + @:analytical() + + if ((q_prim_vf(1)%sf(i, j, 0) < 1.e-10) .and. (model_eqns == 4)) then + !zero density, reassign according to Tait EOS + q_prim_vf(1)%sf(i, j, 0) = & + (((q_prim_vf(E_idx)%sf(i, j, 0) + pi_inf)/(pref + pi_inf))**(1d0/lit_gamma))* & + rhoref*(1d0 - q_prim_vf(alf_idx)%sf(i, j, 0)) + end if + + ! Updating the patch identities bookkeeping variable + if (1d0 - eta < 1d-16) patch_id_fp(i, j, 0) = patch_id end if diff --git a/src/simulation/m_acoustic_src.fpp b/src/simulation/m_acoustic_src.fpp index c0d3963b8..02ee73509 100644 --- a/src/simulation/m_acoustic_src.fpp +++ b/src/simulation/m_acoustic_src.fpp @@ -17,6 +17,8 @@ module m_acoustic_src use m_variables_conversion !< State variables type conversion procedures use m_helper_basic !< Functions to compare floating point numbers + + use m_constants !< Definitions of the constants ! ========================================================================== implicit none private; public :: s_initialize_acoustic_src, s_precalculate_acoustic_spatial_sources, s_acoustic_src_calculations diff --git a/src/simulation/m_boundary_conditions.fpp b/src/simulation/m_boundary_conditions.fpp index ca77632dc..266be8ed0 100644 --- a/src/simulation/m_boundary_conditions.fpp +++ b/src/simulation/m_boundary_conditions.fpp @@ -19,7 +19,7 @@ module m_boundary_conditions implicit none private; - public :: s_populate_primitive_variables_buffers, & + public :: s_populate_variables_buffers, & s_populate_capillary_buffers contains @@ -27,8 +27,7 @@ contains !> The purpose of this procedure is to populate the buffers !! of the primitive variables, depending on the selected !! boundary conditions. - !! @param q_prim_vf Primitive variable - subroutine s_populate_primitive_variables_buffers(q_prim_vf, pb, mv) + subroutine s_populate_variables_buffers(q_prim_vf, pb, mv) type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: pb, mv @@ -213,7 +212,7 @@ contains ! END: Population of Buffers in z-direction ======================== - end subroutine s_populate_primitive_variables_buffers + end subroutine s_populate_variables_buffers subroutine s_ghost_cell_extrapolation(q_prim_vf, pb, mv, bc_dir, bc_loc) diff --git a/src/simulation/m_chemistry.fpp b/src/simulation/m_chemistry.fpp new file mode 100644 index 000000000..d1d9b67d6 --- /dev/null +++ b/src/simulation/m_chemistry.fpp @@ -0,0 +1,188 @@ +!!> +!! @file m_chemistry.f90 +!! @brief Contains module m_chemistry +!! @author Henry Le Berre + +#:include 'macros.fpp' +#:include 'case.fpp' + +module m_chemistry + + use ieee_arithmetic + + use m_mpi_proxy + use m_thermochem + use m_global_parameters + use m_finite_differences + + implicit none + + type(int_bounds_info), private :: ix, iy, iz + type(scalar_field), private :: grads(1:3) + + !$acc declare create(ix, iy, iz) + !$acc declare create(grads) + +contains + + subroutine s_initialize_chemistry_module + + integer :: i + + ix%beg = -buff_size + if (n > 0) then; iy%beg = -buff_size; else; iy%beg = 0; end if + if (p > 0) then; iz%beg = -buff_size; else; iz%beg = 0; end if + + ix%end = m - ix%beg; iy%end = n - iy%beg; iz%end = p - iz%beg + + !$acc update device(ix, iy, iz) + + do i = 1, 3 + @:ALLOCATE(grads(i)%sf(ix%beg:ix%end, iy%beg:iy%end, iz%beg:iz%end)) + end do + + !$acc kernels + do i = 1, num_dims + grads(i)%sf(:, :, :) = 0.0d0 + end do + !$acc end kernels + + end subroutine s_initialize_chemistry_module + + subroutine s_finalize_chemistry_module + + deallocate (grads(1)%sf, grads(2)%sf, grads(3)%sf) + + end subroutine s_finalize_chemistry_module + + #:for NORM_DIR, XYZ in [(1, 'x'), (2, 'y'), (3, 'z')] + subroutine s_compute_chemistry_rhs_${XYZ}$ (flux_n, rhs_vf, flux_src_vf, q_prim_vf) + + type(vector_field), dimension(:), intent(IN) :: flux_n + type(scalar_field), dimension(sys_size), intent(INOUT) :: rhs_vf, flux_src_vf, q_prim_vf + type(int_bounds_info) :: ix, iy, iz + + integer :: x, y, z + integer :: eqn + + integer, parameter :: mx = ${1 if NORM_DIR == 1 else 0}$ + integer, parameter :: my = ${1 if NORM_DIR == 2 else 0}$ + integer, parameter :: mz = ${1 if NORM_DIR == 3 else 0}$ + + !$acc parallel loop collapse(4) present(rhs_vf, flux_n) + do x = 0, m + do y = 0, n + do z = 0, p + + do eqn = chemxb, chemxe + + ! \nabla \cdot (F) + rhs_vf(eqn)%sf(x, y, z) = rhs_vf(eqn)%sf(x, y, z) + & + (flux_n(${NORM_DIR}$)%vf(eqn)%sf(x - mx, y - my, z - mz) - & + flux_n(${NORM_DIR}$)%vf(eqn)%sf(x, y, z))/d${XYZ}$ (${XYZ}$) + + end do + + end do + end do + end do + + end subroutine s_compute_chemistry_rhs_${XYZ}$ + #:endfor + + subroutine s_compute_chemistry_reaction_flux(rhs_vf, q_cons_qp, q_prim_qp) + + type(scalar_field), dimension(sys_size), intent(INOUT) :: rhs_vf, q_cons_qp, q_prim_qp + + integer :: i + + integer :: x, y, z + integer :: eqn + + real(kind(0d0)) :: T + integer :: o + real(kind(0d0)) :: dyn_pres + real(kind(0d0)) :: E + + real(kind(0d0)) :: rho + real(kind(1.d0)), dimension(num_species) :: Ys + real(kind(1.d0)), dimension(num_species) :: omega + real(kind(0d0)), dimension(num_species) :: enthalpies + real(kind(0d0)) :: cp_mix + + #:if chemistry + + !$acc parallel loop collapse(4) private(rho) + do x = 0, m + do y = 0, n + do z = 0, p + + ! Maybe use q_prim_vf instead? + rho = 0d0 + do eqn = chemxb, chemxe + rho = rho + q_cons_qp(eqn)%sf(x, y, z) + end do + + do eqn = chemxb, chemxe + Ys(eqn - chemxb + 1) = q_cons_qp(eqn)%sf(x, y, z)/rho + end do + + dyn_pres = 0d0 + + do i = momxb, momxe + dyn_pres = dyn_pres + rho*q_cons_qp(i)%sf(x, y, z)* & + q_cons_qp(i)%sf(x, y, z)/2d0 + end do + + call get_temperature(.true., q_cons_qp(E_idx)%sf(x, y, z) - dyn_pres, & + & q_prim_qp(tempxb)%sf(x, y, z), Ys, T) + + call get_net_production_rates(rho, T, Ys, omega) + + q_cons_qp(tempxb)%sf(x, y, z) = T + q_prim_qp(tempxb)%sf(x, y, z) = T + + !print*, x, y, z, T, rho, Ys, omega, q_cons_qp(E_idx)%sf(x, y, z), dyn_pres + + do eqn = chemxb, chemxe + + rhs_vf(eqn)%sf(x, y, z) = rhs_vf(eqn)%sf(x, y, z) + & + mol_weights(eqn - chemxb + 1)*omega(eqn - chemxb + 1) + + end do + + end do + end do + end do + + #:else + + @:ASSERT(.false., "Chemistry is not enabled") + + #:endif + + end subroutine s_compute_chemistry_reaction_flux + + subroutine s_chemistry_normalize_cons(q_cons_qp) + + type(scalar_field), dimension(sys_size), intent(INOUT) :: q_cons_qp + + integer :: x, y, z + integer :: eqn + + !$acc parallel loop collapse(4) + do x = ix%beg, ix%end + do y = iy%beg, iy%end + do z = iz%beg, iz%end + + do eqn = chemxb, chemxe + q_cons_qp(eqn)%sf(x, y, z) = max(0d0, q_cons_qp(eqn)%sf(x, y, z)) + end do + + end do + end do + end do + + end subroutine s_chemistry_normalize_cons + +end module m_chemistry diff --git a/src/simulation/m_compute_cbc.fpp b/src/simulation/m_compute_cbc.fpp index c8386c81a..c4b369945 100644 --- a/src/simulation/m_compute_cbc.fpp +++ b/src/simulation/m_compute_cbc.fpp @@ -101,7 +101,7 @@ contains !! any reflections caused by outgoing waves. subroutine s_compute_nonreflecting_subsonic_inflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds) #ifdef CRAY_ACC_WAR - !DIR$ INLINEALWAYS ss_compute_nonreflecting_subsonic_inflow_L + !DIR$ INLINEALWAYS s_compute_nonreflecting_subsonic_inflow_L #else !$acc routine seq #endif diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index d210390dd..b5d33da11 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -348,6 +348,11 @@ contains call s_mpi_abort('ICFL is greater than 1.0. Exiting ...') end if + do i = chemxb, chemxe + !@:ASSERT(all(q_prim_vf(i)%sf(:,:,:) >= -1d0), "bad conc") + !@:ASSERT(all(q_prim_vf(i)%sf(:,:,:) <= 2d0), "bad conc") + end do + if (vcfl_max_glb /= vcfl_max_glb) then call s_mpi_abort('VCFL is NaN. Exiting ...') elseif (vcfl_max_glb > 1d0) then @@ -524,10 +529,8 @@ contains open (2, FILE=trim(file_path)) do j = 0, m - if (((i >= cont_idx%beg) .and. (i <= cont_idx%end)) & - .or. & - ((i >= adv_idx%beg) .and. (i <= adv_idx%end)) & - ) then + ! todo: revisit change here + if (((i >= adv_idx%beg) .and. (i <= adv_idx%end))) then write (2, FMT) x_cb(j), q_cons_vf(i)%sf(j, 0, 0) else write (2, FMT) x_cb(j), q_prim_vf(i)%sf(j, 0, 0) @@ -633,6 +636,8 @@ contains if (((i >= cont_idx%beg) .and. (i <= cont_idx%end)) & .or. & ((i >= adv_idx%beg) .and. (i <= adv_idx%end)) & + .or. & + ((i >= chemxb) .and. (i <= chemxe)) & ) then write (2, FMT) x_cb(j), y_cb(k), q_cons_vf(i)%sf(j, k, 0) else @@ -714,6 +719,8 @@ contains if (((i >= cont_idx%beg) .and. (i <= cont_idx%end)) & .or. & ((i >= adv_idx%beg) .and. (i <= adv_idx%end)) & + .or. & + ((i >= chemxb) .and. (i <= chemxe)) & ) then write (2, FMT) x_cb(j), y_cb(k), z_cb(l), q_cons_vf(i)%sf(j, k, l) else @@ -955,7 +962,7 @@ contains real(kind(0d0)) :: G real(kind(0d0)) :: dyn_p - integer :: i, j, k, l, s, q !< Generic loop iterator + integer :: i, j, k, l, s, q, d !< Generic loop iterator real(kind(0d0)) :: nondim_time !< Non-dimensional time @@ -969,6 +976,8 @@ contains real(kind(0d0)) :: rad, thickness !< For integral quantities logical :: trigger !< For integral quantities + real(kind(0d0)) :: rhoYks(1:num_species) + ! Non-dimensional time calculation if (time_stepper == 23) then nondim_time = mytime @@ -1020,6 +1029,12 @@ contains k = 0 l = 0 + if (chemistry) then + do d = 1, num_species + rhoYks(d) = q_cons_vf(chemxb + d - 1)%sf(j - 2, k, l) + end do + end if + ! Computing/Sharing necessary state variables if (hypoelasticity) then call s_convert_to_mixture_variables(q_cons_vf, j - 2, k, l, & @@ -1039,14 +1054,14 @@ contains call s_compute_pressure( & q_cons_vf(1)%sf(j - 2, k, l), & q_cons_vf(alf_idx)%sf(j - 2, k, l), & - dyn_p, pi_inf, gamma, rho, qv, pres, & + dyn_p, pi_inf, gamma, rho, qv, rhoYks(:), pres, & q_cons_vf(stress_idx%beg)%sf(j - 2, k, l), & q_cons_vf(mom_idx%beg)%sf(j - 2, k, l), G) else call s_compute_pressure( & q_cons_vf(1)%sf(j - 2, k, l), & q_cons_vf(alf_idx)%sf(j - 2, k, l), & - dyn_p, pi_inf, gamma, rho, qv, pres) + dyn_p, pi_inf, gamma, rho, qv, rhoYks(:), pres) end if if (model_eqns == 4) then @@ -1109,6 +1124,12 @@ contains accel = accel_mag(j - 2, k, l) end if elseif (p == 0) then ! 2D simulation + if (chemistry) then + do d = 1, num_species + rhoYks(d) = q_cons_vf(chemxb + d - 1)%sf(j - 2, k - 2, l) + end do + end if + if ((probe(i)%x >= x_cb(-1)) .and. (probe(i)%x <= x_cb(m))) then if ((probe(i)%y >= y_cb(-1)) .and. (probe(i)%y <= y_cb(n))) then do s = -1, m @@ -1139,13 +1160,16 @@ contains call s_compute_pressure( & q_cons_vf(1)%sf(j - 2, k - 2, l), & q_cons_vf(alf_idx)%sf(j - 2, k - 2, l), & - dyn_p, pi_inf, gamma, rho, qv, pres, & + dyn_p, pi_inf, gamma, rho, qv, & + rhoYks, & + pres, & q_cons_vf(stress_idx%beg)%sf(j - 2, k - 2, l), & q_cons_vf(mom_idx%beg)%sf(j - 2, k - 2, l), G) else call s_compute_pressure(q_cons_vf(E_idx)%sf(j - 2, k - 2, l), & q_cons_vf(alf_idx)%sf(j - 2, k - 2, l), & - dyn_p, pi_inf, gamma, rho, qv, pres) + dyn_p, pi_inf, gamma, rho, qv, & + rhoYks, pres) end if if (model_eqns == 4) then @@ -1218,17 +1242,25 @@ contains dyn_p = 0.5d0*rho*dot_product(vel, vel) + if (chemistry) then + do d = 1, num_species + rhoYks(d) = q_cons_vf(chemxb + d - 1)%sf(j - 2, k - 2, l - 2) + end do + end if + if (hypoelasticity) then call s_compute_pressure( & q_cons_vf(1)%sf(j - 2, k - 2, l - 2), & q_cons_vf(alf_idx)%sf(j - 2, k - 2, l - 2), & - dyn_p, pi_inf, gamma, rho, qv, pres, & + dyn_p, pi_inf, gamma, rho, qv, & + rhoYks, pres, & q_cons_vf(stress_idx%beg)%sf(j - 2, k - 2, l - 2), & q_cons_vf(mom_idx%beg)%sf(j - 2, k - 2, l - 2), G) else call s_compute_pressure(q_cons_vf(E_idx)%sf(j - 2, k - 2, l - 2), & q_cons_vf(alf_idx)%sf(j - 2, k - 2, l - 2), & - dyn_p, pi_inf, gamma, rho, qv, pres) + dyn_p, pi_inf, gamma, rho, qv, & + rhoYks, pres) end if ! Compute mixture sound speed diff --git a/src/simulation/m_derived_variables.f90 b/src/simulation/m_derived_variables.f90 index 4cbe11ab3..1c4e838a5 100644 --- a/src/simulation/m_derived_variables.f90 +++ b/src/simulation/m_derived_variables.f90 @@ -21,6 +21,9 @@ module m_derived_variables use m_time_steppers !< Time-stepping algorithms use m_helper + + use m_finite_differences + ! ========================================================================== implicit none diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index 8797a8f45..fda004d79 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -157,6 +157,7 @@ module m_global_parameters logical :: null_weights !< Null undesired WENO weights logical :: mixture_err !< Mixture properties correction logical :: hypoelasticity !< hypoelasticity modeling + logical, parameter :: chemistry = .${chemistry}$. !< Chemistry modeling logical :: cu_tensor logical :: bodyForces @@ -231,6 +232,8 @@ module m_global_parameters integer :: pi_inf_idx !< Index of liquid stiffness func. eqn. type(int_bounds_info) :: stress_idx !< Indexes of first and last shear stress eqns. integer :: c_idx ! Index of the color function + type(int_bounds_info) :: chemistry_idx !< Indexes of first & last concentration eqns. + type(int_bounds_info) :: temperature_idx !< Indexes of first & last temperature eqns. !> @} !$acc declare create(bub_idx) @@ -283,7 +286,7 @@ module m_global_parameters integer :: startx, starty, startz - !$acc declare create(sys_size, buff_size, startx, starty, startz, E_idx, gamma_idx, pi_inf_idx, alf_idx, n_idx, stress_idx) + !$acc declare create(sys_size, buff_size, startx, starty, startz, E_idx, gamma_idx, pi_inf_idx, alf_idx, n_idx, stress_idx, chemistry_idx) ! END: Simulation Algorithm Parameters ===================================== @@ -405,6 +408,9 @@ module m_global_parameters #endif !> @} + type(chemistry_parameters) :: chem_params + !$acc declare create(chem_params) + !> @name Physical bubble parameters (see Ando 2010, Preston 2007) !> @{ real(kind(0d0)) :: R_n, R_v, phi_vn, phi_nv, Pe_c, Tw, pv, M_n, M_v @@ -445,7 +451,10 @@ module m_global_parameters integer :: intxb, intxe integer :: bubxb, bubxe integer :: strxb, strxe -!$acc declare create(momxb, momxe, advxb, advxe, contxb, contxe, intxb, intxe, bubxb, bubxe, strxb, strxe) + integer :: chemxb, chemxe + integer :: tempxb, tempxe + +!$acc declare create(momxb, momxe, advxb, advxe, contxb, contxe, intxb, intxe, bubxb, bubxe, strxb, strxe, chemxb, chemxe, tempxb, tempxe) #ifdef CRAY_ACC_WAR @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:), gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps) @@ -544,6 +553,10 @@ contains teno = .false. #:endif + chem_params%advection = .false. + chem_params%diffusion = .false. + chem_params%reactions = .false. + bc_x%beg = dflt_int; bc_x%end = dflt_int bc_y%beg = dflt_int; bc_y%end = dflt_int bc_z%beg = dflt_int; bc_z%end = dflt_int @@ -1054,6 +1067,16 @@ contains grid_geometry = 3 end if + if (chemistry) then + chemistry_idx%beg = sys_size + 1 + chemistry_idx%end = sys_size + num_species + sys_size = chemistry_idx%end + + temperature_idx%beg = sys_size + 1 + temperature_idx%end = sys_size + 1 + sys_size = temperature_idx%end + end if + momxb = mom_idx%beg momxe = mom_idx%end advxb = adv_idx%beg @@ -1066,8 +1089,13 @@ contains strxe = stress_idx%end intxb = internalEnergies_idx%beg intxe = internalEnergies_idx%end + chemxb = chemistry_idx%beg + chemxe = chemistry_idx%end + tempxb = temperature_idx%beg + tempxe = temperature_idx%end - !$acc update device(momxb, momxe, advxb, advxe, contxb, contxe, bubxb, bubxe, intxb, intxe, sys_size, buff_size, E_idx, alf_idx, n_idx, adv_n, adap_dt, pi_fac, strxb, strxe) + !$acc update device(momxb, momxe, advxb, advxe, contxb, contxe, bubxb, bubxe, intxb, intxe, sys_size, buff_size, E_idx, alf_idx, n_idx, adv_n, adap_dt, pi_fac, strxb, strxe, chemxb, chemxe, tempxb, tempxe) + !$acc update device(chemistry_idx) !$acc update device(cfl_target, m, n, p) !$acc update device(alt_soundspeed, acoustic_source, num_source) diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp index b395fc624..b3d7d53a2 100644 --- a/src/simulation/m_mpi_proxy.fpp +++ b/src/simulation/m_mpi_proxy.fpp @@ -196,8 +196,8 @@ contains call MPI_BCAST(${VAR}$, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) #:endfor - #:for VAR in [ 'run_time_info','cyl_coord', 'mpp_lim', & - & 'mp_weno', 'rdma_mpi', 'weno_flat', 'riemann_flat', & + #:for VAR in [ 'run_time_info','cyl_coord', 'mpp_lim', & + & 'mp_weno', 'rdma_mpi', 'weno_flat', 'riemann_flat', & & 'weno_Re_flux', 'alt_soundspeed', 'null_weights', 'mixture_err', & & 'parallel_io', 'hypoelasticity', 'bubbles', 'polytropic', & & 'polydisperse', 'qbmm', 'acoustic_source', 'probe_wrt', 'integral_wrt', & @@ -207,6 +207,12 @@ contains call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) #:endfor + if (chemistry) then + #:for VAR in [ 'advection', 'diffusion', 'reactions' ] + call MPI_BCAST(chem_params%${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) + #:endfor + end if + #:for VAR in [ 'dt','weno_eps','teno_CT','pref','rhoref','R0ref','Web','Ca', 'sigma', & & 'Re_inv', 'poly_sigma', 'palpha_eps', 'ptgalpha_eps', 'pi_fac', & & 'bc_x%vb1','bc_x%vb2','bc_x%vb3','bc_x%ve1','bc_x%ve2','bc_x%ve2', & diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp index 04693d8fc..2c8b99b5e 100644 --- a/src/simulation/m_rhs.fpp +++ b/src/simulation/m_rhs.fpp @@ -2,6 +2,7 @@ !! @file m_rhs.f90 !! @brief Contains module m_rhs +#:include 'case.fpp' #:include 'macros.fpp' !> @brief The module contains the subroutines used to calculate the right- @@ -54,6 +55,8 @@ module m_rhs use m_surface_tension use m_body_forces + + use m_chemistry ! ========================================================================== implicit none @@ -600,6 +603,16 @@ contains & iz%beg:iz%end)) end do end if + + if (chemistry) then + do l = chemxb, chemxe + @:ALLOCATE(flux_src_n(i)%vf(l)%sf( & + & ix%beg:ix%end, & + & iy%beg:iy%end, & + & iz%beg:iz%end)) + end do + end if + else do l = 1, sys_size @:ALLOCATE(flux_gsrc_n(i)%vf(l)%sf( & @@ -728,7 +741,7 @@ contains real(kind(0d0)) :: start, finish real(kind(0d0)) :: s2, const_sos, s1 - integer :: i, j, k, l, q, ii, id !< Generic loop iterators + integer :: i, c, j, k, l, q, ii, id !< Generic loop iterators integer :: term_index call nvtxStartRange("Compute_RHS") @@ -755,6 +768,13 @@ contains end do end do + if (chemistry) then + !$acc parallel loop default(present) + do i = chemxb, chemxe + rhs_vf(i)%sf(:, :, :) = 0d0 + end do + end if + ! ================================================================== ! Converting Conservative to Primitive Variables ================== @@ -788,7 +808,8 @@ contains call nvtxEndRange call nvtxStartRange("RHS-MPI") - call s_populate_primitive_variables_buffers(q_prim_qp%vf, pb, mv) + call s_populate_variables_buffers(q_prim_qp%vf, pb, mv) + call nvtxEndRange if (cfl_dt) then @@ -975,6 +996,26 @@ contains call nvtxEndRange ! END: Additional physics and source terms ========================= + #:if chemistry + if (chem_params%advection) then + call nvtxStartRange("RHS_Chem_Advection") + + #:for NORM_DIR, XYZ in [(1, 'x'), (2, 'y'), (3, 'z')] + + if (id == ${NORM_DIR}$) then + call s_compute_chemistry_rhs_${XYZ}$ ( & + flux_n, & + rhs_vf, & + flux_src_n(${NORM_DIR}$)%vf, & + q_prim_vf) + end if + + #:endfor + + call nvtxEndRange + end if + #:endif + end do ! END: Dimensional Splitting Loop ================================= @@ -1010,6 +1051,15 @@ contains t_step, & rhs_vf) call nvtxEndRange + + #:if chemistry + if (chem_params%reactions) then + call nvtxStartRange("RHS_Chem_Reactions") + call s_compute_chemistry_reaction_flux(rhs_vf, q_cons_vf, q_prim_qp%vf) + call nvtxEndRange + end if + #:endif + ! END: Additional pphysics and source terms ============================ if (run_time_info .or. probe_wrt .or. ib) then @@ -2384,3 +2434,4 @@ contains end subroutine s_finalize_rhs_module end module m_rhs + diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index 929476b5f..0210882be 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -18,6 +18,7 @@ !! 2) Harten-Lax-van Leer-Contact (HLLC) !! 3) Exact +#:include 'case.fpp' #:include 'macros.fpp' #:include 'inline_riemann.fpp' @@ -35,6 +36,8 @@ module m_riemann_solvers use m_bubbles !< To get the bubble wall pressure function use m_surface_tension !< To get the capilary fluxes + + use m_chemistry ! ========================================================================== implicit none @@ -354,7 +357,7 @@ contains if (norm_dir == ${NORM_DIR}$) then !$acc parallel loop collapse(3) gang vector default(present) private(alpha_rho_L, alpha_rho_R, vel_L, vel_R, alpha_L, alpha_R, vel_avg, tau_e_L, tau_e_R, G_L, G_R, Re_L, Re_R, & - !$acc rho_avg, h_avg, gamma_avg, s_L, s_R, s_S) + !$acc rho_avg, h_avg, gamma_avg, s_L, s_R, s_S, Y_L, Y_R) do l = is3%beg, is3%end do k = is2%beg, is2%end do j = is1%beg, is1%end @@ -733,6 +736,20 @@ contains flux_rs${XYZ}$_vf(j, k, l, contxe) = 0d0 end if end if + + if (chemistry .and. chem_params%advection) then + !$acc loop seq + do i = chemxb, chemxe + Y_L = qL_prim_rs${XYZ}$_vf(j, k, l, i) + Y_R = qR_prim_rs${XYZ}$_vf(j + 1, k, l, i) + + flux_rs${XYZ}$_vf(j, k, l, i) = (s_M*Y_R*rho_R*vel_R(dir_idx(norm_dir)) & + - s_P*Y_L*rho_L*vel_L(dir_idx(norm_dir)) & + + s_M*s_P*(Y_L*rho_L - Y_R*rho_R)) & + /(s_M - s_P) + flux_src_rs${XYZ}$_vf(j, k, l, i) = 0d0 + end do + end if end do end do end do diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 5414bed23..bee897b58 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -38,6 +38,8 @@ module m_start_up use m_rhs !< Right-hand-side (RHS) evaluation procedures + use m_chemistry !< Chemistry module + use m_data_output !< Run-time info & solution data output procedures use m_time_steppers !< Time-stepping algorithms @@ -145,7 +147,7 @@ contains alt_soundspeed, mixture_err, weno_Re_flux, & null_weights, precision, parallel_io, cyl_coord, & rhoref, pref, bubbles, bubble_model, & - R0ref, & + R0ref, chem_params, & #:if not MFC_CASE_OPTIMIZATION nb, mapped_weno, wenoz, teno, weno_order, num_fluids, & #:endif @@ -1048,7 +1050,9 @@ contains real(kind(0d0)), dimension(2) :: Re real(kind(0d0)) :: pres - integer :: i, j, k, l + integer :: i, j, k, l, c + + real(kind(0d0)), dimension(num_species) :: rhoYks do j = 0, m do k = 0, n @@ -1062,8 +1066,14 @@ contains /max(rho, sgm_eps) end do + if (chemistry) then + do c = 1, num_species + rhoYks(c) = v_vf(chemxb + c - 1)%sf(j, k, l) + end do + end if + call s_compute_pressure(v_vf(E_idx)%sf(j, k, l), 0d0, & - dyn_pres, pi_inf, gamma, rho, qv, pres) + dyn_pres, pi_inf, gamma, rho, qv, rhoYks, pres) do i = 1, num_fluids v_vf(i + internalEnergies_idx%beg - 1)%sf(j, k, l) = v_vf(i + adv_idx%beg - 1)%sf(j, k, l)* & @@ -1151,6 +1161,11 @@ contains end if if (relax) call s_infinite_relaxation_k(q_cons_ts(1)%vf) + + if (chemistry) then + call s_chemistry_normalize_cons(q_cons_ts(1)%vf) + end if + ! Time-stepping loop controls t_step = t_step + 1 @@ -1325,6 +1340,8 @@ contains if (hypoelasticity) call s_initialize_hypoelastic_module() if (relax) call s_initialize_phasechange_module() + if (chemistry) call s_initialize_chemistry_module() + call s_initialize_data_output_module() call s_initialize_derived_variables_module() call s_initialize_time_steppers_module() diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index a4ff54a58..865f3f364 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -37,6 +37,8 @@ module m_time_steppers use m_nvtx + use m_thermochem + use m_body_forces ! ========================================================================== @@ -213,6 +215,20 @@ contains @:ACC_SETUP_SFs(q_prim_vf(c_idx)) end if + if (chemistry) then + do i = chemxb, chemxe + @:ALLOCATE(q_prim_vf(i)%sf(ix_t%beg:ix_t%end, & + iy_t%beg:iy_t%end, & + iz_t%beg:iz_t%end)) + @:ACC_SETUP_SFs(q_prim_vf(i)) + end do + + @:ALLOCATE(q_prim_vf(tempxb)%sf(ix_t%beg:ix_t%end, & + iy_t%beg:iy_t%end, & + iz_t%beg:iz_t%end)) + @:ACC_SETUP_SFs(q_prim_vf(tempxb)) + end if + @:ALLOCATE_GLOBAL(pb_ts(1:2)) !Initialize bubble variables pb and mv at all quadrature nodes for all R0 bins if (qbmm .and. (.not. polytropic)) then @@ -306,6 +322,10 @@ contains integer :: i, j, k, l, q!< Generic loop iterator real(kind(0d0)) :: nR3bar + real(kind(0d0)) :: e_mix + + real(kind(0d0)) :: T + real(kind(0d0)), dimension(num_species) :: Ys ! Stage 1 of 1 ===================================================== diff --git a/src/simulation/m_viscous.fpp b/src/simulation/m_viscous.fpp index 073c7d278..c8c83f6b4 100644 --- a/src/simulation/m_viscous.fpp +++ b/src/simulation/m_viscous.fpp @@ -14,6 +14,8 @@ module m_viscous use m_weno use m_helper + + use m_finite_differences ! ========================================================================== private; public s_get_viscous, & @@ -522,7 +524,7 @@ contains end if end subroutine s_compute_viscous_stress_tensor -!> Computes viscous terms + !> Computes viscous terms !! @param q_cons_vf Cell-averaged conservative variables !! @param q_prim_vf Cell-averaged primitive variables !! @param rhs_vf Cell-averaged RHS variables @@ -964,6 +966,7 @@ contains end if else + do i = iv%beg, iv%end call s_compute_fd_gradient(q_prim_qp%vf(i), & dq_prim_dx_qp(1)%vf(i), & @@ -1015,7 +1018,6 @@ contains if (n > 0) then if (p > 0) then - call s_weno(v_vf(iv%beg:iv%end), & vL_x(:, :, :, iv%beg:iv%end), vL_y(:, :, :, iv%beg:iv%end), vL_z(:, :, :, iv%beg:iv%end), vR_x(:, :, :, iv%beg:iv%end), vR_y(:, :, :, iv%beg:iv%end), vR_z(:, :, :, iv%beg:iv%end), & norm_dir, weno_dir, & @@ -1027,7 +1029,6 @@ contains is1_viscous, is2_viscous, is3_viscous) end if else - call s_weno(v_vf(iv%beg:iv%end), & vL_x(:, :, :, iv%beg:iv%end), vL_y(:, :, :, :), vL_z(:, :, :, :), vR_x(:, :, :, iv%beg:iv%end), vR_y(:, :, :, :), vR_z(:, :, :, :), & norm_dir, weno_dir, & diff --git a/toolchain/mfc/build.py b/toolchain/mfc/build.py index bdcaf32ef..2e44fe9ec 100644 --- a/toolchain/mfc/build.py +++ b/toolchain/mfc/build.py @@ -40,6 +40,9 @@ def get_slug(self, case: input.MFCInputFile) -> str: m.update(CFG().make_slug().encode()) m.update(case.get_fpp(self, False).encode()) + if case.params.get('chemistry', 'F') == 'T': + m.update(case.get_cantera_solution().name.encode()) + return m.hexdigest()[:10] # Get path to directory that will store the build files diff --git a/toolchain/mfc/case.py b/toolchain/mfc/case.py index b62282bb2..12209acbd 100644 --- a/toolchain/mfc/case.py +++ b/toolchain/mfc/case.py @@ -46,7 +46,7 @@ def get_inp(self, _target) -> str: # Create Fortran-style input file content string dict_str = "" for key, val in self.params.items(): - if key in MASTER_KEYS: + if key in MASTER_KEYS and key not in case_dicts.IGNORE: if self.__is_ic_analytical(key, val): dict_str += f"{key} = 0d0\n" ignored.append(key) @@ -211,6 +211,11 @@ def __get_sim_fpp(self, print: bool) -> str: """ def get_fpp(self, target, print = True) -> str: + def _prepend() -> str: + return f"""\ +#:set chemistry = {self.params.get("chemistry", 'F') == 'T'} +""" + def _default(_) -> str: return "! This file is purposefully empty." @@ -219,7 +224,7 @@ def _default(_) -> str: "simulation" : self.__get_sim_fpp, }.get(build.get_target(target).name, _default)(print) - return result + return _prepend() + result def __getitem__(self, key: str) -> str: return self.params[key] diff --git a/toolchain/mfc/run/case_dicts.py b/toolchain/mfc/run/case_dicts.py index 95b2a5bcc..db4e19fdd 100644 --- a/toolchain/mfc/run/case_dicts.py +++ b/toolchain/mfc/run/case_dicts.py @@ -49,6 +49,8 @@ def analytic(self): 'adv_n': ParamType.LOG, 'cfl_adap_dt': ParamType.LOG, 'cfl_const_dt': ParamType.LOG, + 'chemistry': ParamType.LOG, + 'cantera_file': ParamType.STR, } PRE_PROCESS = COMMON.copy() @@ -125,9 +127,8 @@ def analytic(self): PRE_PROCESS[f"patch_icpp({p_id})%{real_attr}"] = ParamType.REAL PRE_PROCESS[f"patch_icpp({p_id})%pres"] = ParamType.REAL.analytic() - # (cameron): This parameter has since been removed. - # for i in range(100): - # PRE_PROCESS.append(f"patch_icpp({p_id})%Y({i})") + for i in range(100): + PRE_PROCESS[f"patch_icpp({p_id})%Y({i})"] = ParamType.REAL.analytic() PRE_PROCESS[f"patch_icpp({p_id})%model%filepath"] = ParamType.STR @@ -221,9 +222,8 @@ def analytic(self): 'low_Mach': ParamType.INT, }) -# NOTE: Not currently present -# for var in [ 'advection', 'diffusion', 'reactions' ]: -# SIMULATION.append(f'chem_params%{var}') +for var in [ 'advection', 'diffusion', 'reactions' ]: + SIMULATION[f'chem_params%{var}'] = ParamType.LOG for ib_id in range(1, 10+1): for real_attr, ty in [("geometry", ParamType.INT), ("radius", ParamType.REAL), @@ -255,12 +255,6 @@ def analytic(self): for prepend in ["domain%beg", "domain%end"]: SIMULATION[f"{cmp}_{prepend}"] = ParamType.REAL -# NOTE: This is now just "probe_wrt" -# for wrt_id in range(1,10+1): -# for cmp in ["x", "y", "z"]: -# SIMULATION.append(f'probe_wrt({wrt_id})%{cmp}') -# set_type(f'probe_wrt({wrt_id})%{cmp}', ParamType.LOG) - for probe_id in range(1,3+1): for cmp in ["x", "y", "z"]: SIMULATION[f'probe({probe_id})%{cmp}'] = ParamType.REAL @@ -295,7 +289,7 @@ def analytic(self): SIMULATION[f"integral({int_id})%{cmp}max"] = ParamType.REAL -# Removed: 'fourier_modes%beg', 'fourier_modes%end', 'chem_wrt' +# Removed: 'fourier_modes%beg', 'fourier_modes%end', 'chem_wrt_Y' # Feel free to return them if they are needed once more. POST_PROCESS = COMMON.copy() POST_PROCESS.update({ @@ -345,9 +339,9 @@ def analytic(self): for real_attr in ["mom_wrt", "vel_wrt", "flux_wrt", "omega_wrt"]: POST_PROCESS[f'{real_attr}({cmp_id})'] = ParamType.LOG -# NOTE: `chem_wrt` is missing +# NOTE: `chem_wrt_Y` is missing # for cmp_id in range(100): -# POST_PROCESS.append(f'chem_wrt({cmp_id})') +# POST_PROCESS.append(f'chem_wrt_Y({cmp_id})') for fl_id in range(1,10+1): for append, ty in [("schlieren_alpha", ParamType.REAL), @@ -359,6 +353,7 @@ def analytic(self): "cv", "qv", "qvp" ]: POST_PROCESS[f"fluid_pp({fl_id})%{real_attr}"] = ParamType.REAL +IGNORE = ["cantera_file", "chemistry"] ALL = COMMON.copy() ALL.update(PRE_PROCESS) diff --git a/toolchain/mfc/run/input.py b/toolchain/mfc/run/input.py index 27110e1a2..b361e0e2e 100644 --- a/toolchain/mfc/run/input.py +++ b/toolchain/mfc/run/input.py @@ -1,5 +1,8 @@ import os, json, glob, typing, dataclasses +import pyrometheus as pyro +import cantera as ct + from ..printer import cons from .. import common, build from ..state import ARGS @@ -30,6 +33,23 @@ def __save_fpp(self, target, contents: str) -> None: cons.print("Writing a (new) custom case.fpp file.") common.file_write(fpp_path, contents, True) + def get_cantera_solution(self) -> ct.Solution: + cantera_file = self.params["cantera_file"] + + candidates = [ + cantera_file, + os.path.join(self.dirpath, cantera_file), + os.path.join(common.MFC_MECHANISMS_DIR, cantera_file), + ] + + for candidate in candidates: + try: + return ct.Solution(candidate) + except Exception: + continue + + raise common.MFCException(f"Cantera file '{cantera_file}' not found. Searched: {', '.join(candidates)}.") + def generate_fpp(self, target) -> None: if target.isDependency: return @@ -40,6 +60,20 @@ def generate_fpp(self, target) -> None: # Case FPP file self.__save_fpp(target, self.get_fpp(target)) + # Chemistry Rates FPP file + modules_dir = os.path.join(target.get_staging_dirpath(self), "modules", target.name) + common.create_directory(modules_dir) + + if self.params.get("chemistry", 'F') == 'T': + common.file_write( + os.path.join(modules_dir, "m_pyrometheus.f90"), + pyro.codegen.fortran90.gen_thermochem_code( + self.get_cantera_solution(), + module_name="m_pyrometheus" + ), + True + ) + cons.unindent() # Generate case.fpp & [target.name].inp diff --git a/toolchain/requirements.txt b/toolchain/requirements.txt index fae4c35c6..e438efa24 100644 --- a/toolchain/requirements.txt +++ b/toolchain/requirements.txt @@ -1,14 +1,28 @@ -numpy -pandas +# General rich -fypp -mako wheel typos typing PyYAML -pylint argparse -fprettify dataclasses jsonschema + +# Build System +fypp + +# Running Cases +mako + +# Code Health +typos +pylint +fprettify + +# Profiling +numpy +pandas + +# Chemistry +cantera +pyrometheus @ git+https://github.com/henryleberre/pyrometheus.git@multi-lang-part-2