diff --git a/.vscode/launch.json b/.vscode/launch.json index 76f40d7..14e8f1f 100644 --- a/.vscode/launch.json +++ b/.vscode/launch.json @@ -43,5 +43,19 @@ "externalConsole": false, "MIMode": "gdb", }, + { + "name": "(gdb) Sedov 2D", + "type": "cppdbg", + "request": "launch", + "program": "${workspaceFolder}/build/bin/cato.x", + "args": [ + "input.ini" + ], + "stopAtEntry": false, + "cwd": "${workspaceFolder}/tests/integrated/sedov", + "environment": [], + "externalConsole": false, + "MIMode": "gdb", + }, ] } diff --git a/src/lib/fluid/mod_fluid.f90 b/src/lib/fluid/mod_fluid.f90 index be3acb3..4f467be 100644 --- a/src/lib/fluid/mod_fluid.f90 +++ b/src/lib/fluid/mod_fluid.f90 @@ -1,5 +1,5 @@ ! MIT License -! Copyright (c) 2019 Sam Miller +! Copyright (c) 2020 Sam Miller ! Permission is hereby granted, free of charge, to any person obtaining a copy ! of this software and associated documentation files (the "Software"), to deal ! in the Software without restriction, including without limitation the rights @@ -541,7 +541,7 @@ subroutine get_continuity_sensor(self) !< or discontinuous (linear or non-linear) class(fluid_t), intent(inout) :: self - call distinguish(rho=self%rho, u=self%u, v=self%v, p=self%p, continuity_sensor=self%continuous_sensor) + call distinguish(lbounds=lbound(self%rho), rho=self%rho, u=self%u, v=self%v, p=self%p, continuity_sensor=self%continuous_sensor) end subroutine get_continuity_sensor subroutine residual_smoother(self) diff --git a/src/lib/io/contour_writer.f90 b/src/lib/io/contour_writer.f90 index efbb7bf..2a4d6a7 100644 --- a/src/lib/io/contour_writer.f90 +++ b/src/lib/io/contour_writer.f90 @@ -284,8 +284,8 @@ subroutine write_hdf5(self, master, time, iteration) description='Cell j Index', units='dimensionless') deallocate(int_data_buffer) - call self%write_2d_integer_data(data=master%fluid%continuous_sensor, name='/continuity_sensor', & - description='Continuity Sensor [0=continuous, 1=linear discontinuity, 2=non-linear discontinuity', & + call self%write_2d_integer_data(data=master%fluid%continuous_sensor(ilo:ihi, jlo:jhi), name='/continuity_sensor', & + description='Continuity Sensor [0=continuous, 1=linear discontinuity, 2=non-linear discontinuity]', & units='dimensionless') ! Primitive Variables diff --git a/src/lib/io/mod_input.f90 b/src/lib/io/mod_input.f90 index 0fd3cb2..fcd9b76 100644 --- a/src/lib/io/mod_input.f90 +++ b/src/lib/io/mod_input.f90 @@ -287,23 +287,23 @@ subroutine read_from_ini(self, filename) call cfg%get("initial_conditions", "init_pressure", self%init_pressure, 0.0_rk) end if - ! Grid - select case(trim(self%edge_interpolation_scheme)) - case('TVD2', 'TVD3', 'TVD5', 'MLP3', 'MLP5') - required_n_ghost_layers = 2 - call cfg%get("grid", "n_ghost_layers", self%n_ghost_layers, required_n_ghost_layers) - case default - call error_msg(module='mod_input', class="input_t", procedure='read_from_ini', & - message="Unknown edge interpolation scheme, must be one of the following: "// & - "'TVD2', 'TVD3', 'TVD5', 'MLP3', or 'MLP5'", & - file_name=__FILE__, line_number=__LINE__) - end select - - if(self%n_ghost_layers /= required_n_ghost_layers) then - call error_msg(module='mod_input', class="input_t", procedure='read_from_ini', & - message="The number of required ghost cell layers doesn't match the edge interpolation order", & - file_name=__FILE__, line_number=__LINE__) - end if + ! ! Grid + ! select case(trim(self%edge_interpolation_scheme)) + ! case('TVD2', 'TVD3', 'TVD5', 'MLP3', 'MLP5') + ! required_n_ghost_layers = 2 + ! call cfg%get("grid", "n_ghost_layers", self%n_ghost_layers, required_n_ghost_layers) + ! case default + ! call error_msg(module='mod_input', class="input_t", procedure='read_from_ini', & + ! message="Unknown edge interpolation scheme, must be one of the following: "// & + ! "'TVD2', 'TVD3', 'TVD5', 'MLP3', or 'MLP5'", & + ! file_name=__FILE__, line_number=__LINE__) + ! end select + + ! if(self%n_ghost_layers /= required_n_ghost_layers) then + ! call error_msg(module='mod_input', class="input_t", procedure='read_from_ini', & + ! message="The number of required ghost cell layers doesn't match the edge interpolation order", & + ! file_name=__FILE__, line_number=__LINE__) + ! end if call cfg%get("grid", "grid_type", char_buffer, '2d_regular') self%grid_type = trim(char_buffer) diff --git a/src/lib/limiters/e_mlp_distinguisher.fypp b/src/lib/limiters/e_mlp_distinguisher.fypp index 5f226f7..6e70f62 100644 --- a/src/lib/limiters/e_mlp_distinguisher.fypp +++ b/src/lib/limiters/e_mlp_distinguisher.fypp @@ -18,6 +18,12 @@ ! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE ! SOFTWARE. +#ifdef __SIMD_ALIGN_OMP__ +#define __CONT_ALIGN__ aligned(rho, u, v, p, d_bar_rho, d_bar_u, d_bar_v, d_bar_p:__ALIGNBYTES__) +#else +#define __CONT_ALIGN__ +#endif + module mod_distinguisher !< Summary: Provide the procedures to conduct the "distinguisher" step in e-MLP. This process !< scans the domain and finds the regions of linear/non-linear discontinuities as well @@ -39,12 +45,14 @@ module mod_distinguisher contains - subroutine distinguish(rho, u, v, p, continuity_sensor) + subroutine distinguish(lbounds, rho, u, v, p, continuity_sensor) !< Scan the domain for continuous and discontinuous regions - real(rk), dimension(:, :), contiguous, intent(in) :: rho !< (i,j); cell-centered density - real(rk), dimension(:, :), contiguous, intent(in) :: u !< (i,j); cell-centered x-velocity - real(rk), dimension(:, :), contiguous, intent(in) :: v !< (i,j); cell-centered y-velocity - real(rk), dimension(:, :), contiguous, intent(in) :: p !< (i,j); cell-centered pressure + + integer(ik), dimension(2), intent(in) :: lbounds + real(rk), dimension(lbounds(1):, lbounds(2):), contiguous, intent(in) :: rho !< (i,j); cell-centered density + real(rk), dimension(lbounds(1):, lbounds(2):), contiguous, intent(in) :: u !< (i,j); cell-centered x-velocity + real(rk), dimension(lbounds(1):, lbounds(2):), contiguous, intent(in) :: v !< (i,j); cell-centered y-velocity + real(rk), dimension(lbounds(1):, lbounds(2):), contiguous, intent(in) :: p !< (i,j); cell-centered pressure integer(ik), dimension(:, :), allocatable, intent(out) :: continuity_sensor !< (i,j); sensor value to tag wheter it is continuous, linear discontinuous, or non-linear discontinuous @@ -70,12 +78,18 @@ contains jhi = ubound(rho, dim=2) allocate(continuity_sensor(ilo:ihi, jlo:jhi)) + continuity_sensor = CONTINUOUS_REGION allocate(d_bar_rho(ilo:ihi, jlo:jhi)) + !dir$ assume_aligned d_bar_rho: __ALIGNBYTES__ allocate(d_bar_u(ilo:ihi, jlo:jhi)) + !dir$ assume_aligned d_bar_u: __ALIGNBYTES__ allocate(d_bar_v(ilo:ihi, jlo:jhi)) + !dir$ assume_aligned d_bar_v: __ALIGNBYTES__ allocate(d_bar_p(ilo:ihi, jlo:jhi)) + !dir$ assume_aligned d_bar_p: __ALIGNBYTES__ + d_bar_rho = 0.0_rk d_bar_u = 0.0_rk d_bar_v = 0.0_rk @@ -89,18 +103,17 @@ contains if(abs(${F}$(i, j)) > 0.0_rk) then ! Eq 11a d_ij_i = abs(((-one_sixth * ${F}$(i - 2, j) + two_thirds * ${F}$(i - 1, j) + two_thirds * ${F}$(i + 1, j) - one_sixth * ${F}$(i + 2, j)) & - / ${F}$(i,j)) & - - 1.0_rk) + / ${F}$(i,j)) - 1.0_rk) ! Eq 11b d_ij_j = abs(((-one_sixth * ${F}$(i, j - 2) + two_thirds * ${F}$(i, j - 1) + two_thirds * ${F}$(i, j + 1) - one_sixth * ${F}$(i, j + 2)) & - / ${F}$(i,j)) & - - 1.0_rk) + / ${F}$(i,j)) - 1.0_rk) end if - d_bar_${F}$(i, j) = 0.5_rk * (d_ij_i + d_ij_j) ! Eq 11c + d_bar_${F}$(i, j) = 0.5_rk * (d_ij_i + d_ij_j) end do end do + #:endfor ! Assign contiuous sensor based on the approximate values, e.g. d_bar_rho @@ -108,16 +121,18 @@ contains do i = ilo, ihi if(abs(d_bar_rho(i, j)) > EPS) continuity_sensor(i, j) = LINEAR_DISCONT_REGION - if(abs(u(i, j)) < abs(v(i, j))) then - if(abs(d_bar_v(i, j)) > EPS .and. abs(v(i, j)) > 1e-6_rk) continuity_sensor(i, j) = LINEAR_DISCONT_REGION - else ! abs(u(i,j)) >= abs(v(i,j)) + if(abs(u(i, j)) > abs(v(i, j)) .or. abs(u(i, j) - v(i, j)) < epsilon(1.0_rk)) then if(abs(d_bar_u(i, j)) > EPS .and. abs(u(i, j)) > 1e-6_rk) continuity_sensor(i, j) = LINEAR_DISCONT_REGION + + else if(abs(u(i, j)) < abs(v(i, j))) then + if(abs(d_bar_v(i, j)) > EPS .and. abs(v(i, j)) > 1e-6_rk) continuity_sensor(i, j) = LINEAR_DISCONT_REGION end if if(abs(d_bar_p(i, j)) > EPS) continuity_sensor(i, j) = NONLINEAR_DISCONT_REGION end do end do + deallocate(d_bar_rho) deallocate(d_bar_u) deallocate(d_bar_v) diff --git a/src/lib/spatial_reconstruction/muscl_e_mlp.fypp b/src/lib/spatial_reconstruction/muscl_e_mlp.fypp index d44a8fe..2500375 100644 --- a/src/lib/spatial_reconstruction/muscl_e_mlp.fypp +++ b/src/lib/spatial_reconstruction/muscl_e_mlp.fypp @@ -146,7 +146,7 @@ subroutine distinguish_continuous_regions(self, rho, u, v, p, lbounds) real(rk) :: d_ij_i, d_ij_j - real(rk), parameter :: EPS = 0.01_rk + real(rk), parameter :: EPS = 0.005_rk ! Ref [1] Recommends 0.001 for steady problems and 0.01 for other cases... real(rk), parameter :: one_sixth = 1.0_rk / 6.0_rk @@ -231,6 +231,8 @@ subroutine distinguish_continuous_regions(self, rho, u, v, p, lbounds) deallocate(d_bar_v) deallocate(d_bar_p) + self%continuity_sensor = LINEAR_DISCONT_REGION + end subroutine distinguish_continuous_regions subroutine finalize(self) @@ -255,22 +257,22 @@ subroutine interp_muscl_e_mlp_${DIR}$_edge_${ORDER}$_order(self, q, lbounds, edg integer(ik) :: i, j, m, n integer(ik) :: ilo, ihi, jlo, jhi - real(rk) :: delta_i_minus_half !< Delta Phi_{i-1/2, j} = q(i, j) - q(i - 1, j) (or j-1/2 depending on direction) - real(rk) :: delta_i_plus_half !< Delta Phi_{i+1/2, j} = q(i + 1, j) - q(i, j) (or j+1/2 depending on direction) - real(rk) :: delta_i_plus_three_half !< Delta Phi_{i+3/2, j} = q(i + 2, j) - q(i + 1, j) (or j+3/2 depending on direction) + real(rk) :: delta_ij_minus_half !< Delta Phi_{i-1/2, j} = q(i, j) - q(i - 1, j) (or j-1/2 depending on direction) + real(rk) :: delta_ij_plus_half !< Delta Phi_{i+1/2, j} = q(i + 1, j) - q(i, j) (or j+1/2 depending on direction) + real(rk) :: delta_ij_plus_three_half !< Delta Phi_{i+3/2, j} = q(i + 2, j) - q(i + 1, j) (or j+3/2 depending on direction) real(rk) :: r_R !< smoothness parameter real(rk) :: r_L !< smoothness parameter - real(rk) :: tan_theta_i !< the tan(theta), i term - real(rk) :: tan_theta_ip1 !< the tan(theta), i+1 term + real(rk) :: tan_theta_ij !< the tan(theta), i term + real(rk) :: tan_theta_ijp1 !< the tan(theta), i+1 term real(rk) :: phi_ML, phi_MR real(rk), dimension(6) :: phi_MN real(rk) :: min_phi_MN, max_phi_MN real(rk) :: alpha_L, alpha_R real(rk) :: alpha_L_term, alpha_R_term - real(rk) :: beta_L_delta_i_minus_half - real(rk) :: beta_R_delta_i_plus_three_half + real(rk) :: beta_L_delta_ij_minus_half + real(rk) :: beta_R_delta_ij_plus_three_half real(rk), parameter :: f_1_over_3 = 1.0_rk / 3.0_rk !< paramater for 1/3 used to save extra divisions real(rk), parameter :: f_2_over_3 = 2.0_rk / 3.0_rk !< paramater for 2/3 used to save extra divisions @@ -314,11 +316,11 @@ subroutine interp_muscl_e_mlp_${DIR}$_edge_${ORDER}$_order(self, q, lbounds, edg !$omp parallel default(none), & !$omp firstprivate(ilo, ihi, jlo, jhi) & - !$omp private(i, j, m, n, delta_i_minus_half, delta_i_plus_half, delta_i_plus_three_half) & + !$omp private(i, j, m, n, delta_ij_minus_half, delta_ij_plus_half, delta_ij_plus_three_half) & !$omp private(r_L, r_R) & - !$omp private(tan_theta_i, tan_theta_ip1, alpha_L, alpha_R, alpha_L_term, alpha_R_term) & + !$omp private(tan_theta_ij, tan_theta_ijp1, alpha_L, alpha_R, alpha_L_term, alpha_R_term) & !$omp private(phi_ML, phi_MR) & - !$omp private(beta_L_delta_i_minus_half, beta_R_delta_i_plus_three_half, phi_MN, min_phi_MN, max_phi_MN) & + !$omp private(beta_L_delta_ij_minus_half, beta_R_delta_ij_plus_three_half, phi_MN, min_phi_MN, max_phi_MN) & !$omp shared(q, edge_values, self) !$omp do do j = jlo, jhi @@ -329,36 +331,36 @@ subroutine interp_muscl_e_mlp_${DIR}$_edge_${ORDER}$_order(self, q, lbounds, edg ! 3rd order interpolation #:if DIR == 'i' ! Eq 23a in Ref [1] - beta_L_delta_i_minus_half = -f_1_over_3 * q(i - 1, j) - f_1_over_3 * q(i,j) + f_2_over_3 * q(i + 1, j) + beta_L_delta_ij_minus_half = -f_1_over_3 * q(i - 1, j) - f_1_over_3 * q(i,j) + f_2_over_3 * q(i + 1, j) ! Eq 23b in Ref [1] - beta_R_delta_i_plus_three_half = f_1_over_3 * q(i + 2, j) + f_1_over_3 * q(i + 1,j) - f_2_over_3 * q(i, j) + beta_R_delta_ij_plus_three_half = f_1_over_3 * q(i + 2, j) + f_1_over_3 * q(i + 1, j) - f_2_over_3 * q(i, j) #:elif DIR == 'j' ! Eq 23a in Ref [1] - beta_L_delta_i_minus_half = -f_1_over_3 * q(i, j - 1) - f_1_over_3 * q(i,j) + f_2_over_3 * q(i, j + 1) + beta_L_delta_ij_minus_half = -f_1_over_3 * q(i, j - 1) - f_1_over_3 * q(i,j) + f_2_over_3 * q(i, j + 1) ! Eq 23b in Ref [1] - beta_R_delta_i_plus_three_half = f_1_over_3 * q(i, j + 2) + f_1_over_3 * q(i,j + 1) - f_2_over_3 * q(i, j) + beta_R_delta_ij_plus_three_half = f_1_over_3 * q(i, j + 2) + f_1_over_3 * q(i, j + 1) - f_2_over_3 * q(i, j) #:endif #:elif ORDER == 5 ! 5th order interpolation #:if DIR == 'i' ! Eq 24a in Ref [1] - beta_L_delta_i_minus_half = f_2_over_30 * q(i - 2, j) - f_13_over_30 * q(i - 1,j) & + beta_L_delta_ij_minus_half = f_2_over_30 * q(i - 2, j) - f_13_over_30 * q(i - 1,j) & - f_13_over_30 * q(i, j) + f_27_over_30 * q(i + 1, j) - f_3_over_30 * q(i + 2, j) ! Eq 24b in Ref [1] - beta_R_delta_i_plus_three_half = f_2_over_30 * q(i - 2, j) - f_13_over_30 * q(i - 1,j) & - - f_13_over_30 * q(i, j) + f_27_over_30 * q(i + 1, j) - f_3_over_30 * q(i + 2, j) + beta_R_delta_ij_plus_three_half = -f_2_over_30 * q(i + 3, j) + f_13_over_30 * q(i + 2,j) & + + f_13_over_30 * q(i + 1, j) - f_27_over_30 * q(i, j) + f_3_over_30 * q(i - 1, j) #:elif DIR == 'j' ! Eq 24a in Ref [1] - beta_L_delta_i_minus_half = f_2_over_30 * q(i, j - 2) - f_13_over_30 * q(i,j - 1) & + beta_L_delta_ij_minus_half = f_2_over_30 * q(i, j - 2) - f_13_over_30 * q(i, j - 1) & - f_13_over_30 * q(i, j) + f_27_over_30 * q(i, j + 1) - f_3_over_30 * q(i, j + 2) ! Eq 24b in Ref [1] - beta_R_delta_i_plus_three_half = f_2_over_30 * q(i, j - 2) - f_13_over_30 * q(i,j - 1) & - - f_13_over_30 * q(i, j) + f_27_over_30 * q(i, j + 1) - f_3_over_30 * q(i, j + 2) + beta_R_delta_ij_plus_three_half = -f_2_over_30 * q(i, j + 3) + f_13_over_30 * q(i,j + 2) & + + f_13_over_30 * q(i, j + 1) - f_27_over_30 * q(i, j) + f_3_over_30 * q(i, j - 1) #:endif #:endif @@ -368,103 +370,106 @@ subroutine interp_muscl_e_mlp_${DIR}$_edge_${ORDER}$_order(self, q, lbounds, edg ! No limiting function enabled, since it's a continuous region ! Phi L, e.g. left state - edge_values(1, i, j) = q(i, j) + 0.5_rk * beta_L_delta_i_minus_half + edge_values(1, i, j) = q(i, j) + 0.5_rk * beta_L_delta_ij_minus_half #:if DIR == 'i' ! Phi R, e.g. right state - edge_values(2, i, j) = q(i + 1, j) - 0.5_rk * beta_R_delta_i_plus_three_half + edge_values(2, i, j) = q(i + 1, j) - 0.5_rk * beta_R_delta_ij_plus_three_half #:elif DIR == 'j' ! Phi R, e.g. right state - edge_values(2, i, j) = q(i, j + 1) - 0.5_rk * beta_R_delta_i_plus_three_half + edge_values(2, i, j) = q(i, j + 1) - 0.5_rk * beta_R_delta_ij_plus_three_half #:endif case(LINEAR_DISCONT_REGION) ! For a linear discontinuity, use TVD limiting #:if DIR == 'i' - delta_i_minus_half = q(i, j) - q(i - 1, j) - delta_i_plus_half = q(i + 1, j) - q(i, j) - delta_i_plus_three_half = q(i + 2, j) - q(i + 1, j) + delta_ij_minus_half = q(i, j) - q(i - 1, j) + delta_ij_plus_half = q(i + 1, j) - q(i, j) + delta_ij_plus_three_half = q(i + 2, j) - q(i + 1, j) #:elif DIR == 'j' - delta_i_minus_half = q(i, j) - q(i, j - 1) - delta_i_plus_half = q(i, j + 1) - q(i, j) - delta_i_plus_three_half = q(i, j + 2) - q(i, j + 1) + delta_ij_minus_half = q(i, j) - q(i, j - 1) + delta_ij_plus_half = q(i, j + 1) - q(i, j) + delta_ij_plus_three_half = q(i, j + 2) - q(i, j + 1) #:endif - if (abs(delta_i_minus_half) < epsilon(1.0_rk)) delta_i_minus_half = 0.0_rk - if (abs(delta_i_plus_half) < epsilon(1.0_rk)) delta_i_plus_half = 0.0_rk - if (abs(delta_i_plus_three_half) < epsilon(1.0_rk)) delta_i_plus_three_half = 0.0_rk + if (abs(delta_ij_minus_half) < epsilon(1.0_rk)) delta_ij_minus_half = 0.0_rk + if (abs(delta_ij_plus_half) < epsilon(1.0_rk)) delta_ij_plus_half = 0.0_rk + if (abs(delta_ij_plus_three_half) < epsilon(1.0_rk)) delta_ij_plus_three_half = 0.0_rk ! Phi L, e.g. left state - edge_values(1, i, j) = q(i, j) + 0.5_rk * max(0.0_rk, min(2.0_rk * delta_i_plus_half, & - 2.0_rk * delta_i_minus_half, & - beta_L_delta_i_minus_half)) + edge_values(1, i, j) = q(i, j) + 0.5_rk * max(0.0_rk, min(2.0_rk * delta_ij_plus_half, & + 2.0_rk * delta_ij_minus_half, & + beta_L_delta_ij_minus_half)) ! Phi R, e.g. right state #:if DIR == 'i' - edge_values(2, i, j) = q(i + 1, j) - 0.5_rk * max(0.0_rk, min(2.0_rk * delta_i_plus_half, & - 2.0_rk * delta_i_plus_three_half, & - beta_R_delta_i_plus_three_half)) + edge_values(2, i, j) = q(i + 1, j) - 0.5_rk * max(0.0_rk, min(2.0_rk * delta_ij_plus_half, & + 2.0_rk * delta_ij_plus_three_half, & + beta_R_delta_ij_plus_three_half)) #:elif DIR == 'j' - edge_values(2, i, j) = q(i, j + 1) - 0.5_rk * max(0.0_rk, min(2.0_rk * delta_i_plus_half, & - 2.0_rk * delta_i_plus_three_half, & - beta_R_delta_i_plus_three_half)) + edge_values(2, i, j) = q(i, j + 1) - 0.5_rk * max(0.0_rk, min(2.0_rk * delta_ij_plus_half, & + 2.0_rk * delta_ij_plus_three_half, & + beta_R_delta_ij_plus_three_half)) #:endif case(NONLINEAR_DISCONT_REGION) ! For non-linear discontinuity, use MLP limiting #:if DIR == 'i' - delta_i_minus_half = q(i, j) - q(i - 1, j) - delta_i_plus_half = q(i + 1, j) - q(i, j) - delta_i_plus_three_half = q(i + 2, j) - q(i + 1, j) + delta_ij_minus_half = q(i, j) - q(i - 1, j) + delta_ij_plus_half = q(i + 1, j) - q(i, j) + delta_ij_plus_three_half = q(i + 2, j) - q(i + 1, j) #:elif DIR == 'j' - delta_i_minus_half = q(i, j) - q(i, j - 1) - delta_i_plus_half = q(i, j + 1) - q(i, j) - delta_i_plus_three_half = q(i, j + 2) - q(i, j + 1) + delta_ij_minus_half = q(i, j) - q(i, j - 1) + delta_ij_plus_half = q(i, j + 1) - q(i, j) + delta_ij_plus_three_half = q(i, j + 2) - q(i, j + 1) #:endif - if (abs(delta_i_minus_half) < epsilon(1.0_rk)) delta_i_minus_half = 0.0_rk - if (abs(delta_i_plus_half) < epsilon(1.0_rk)) delta_i_plus_half = 0.0_rk - if (abs(delta_i_plus_three_half) < epsilon(1.0_rk)) delta_i_plus_three_half = 0.0_rk + if (abs(delta_ij_minus_half) < epsilon(1.0_rk)) delta_ij_minus_half = 0.0_rk + if (abs(delta_ij_plus_half) < epsilon(1.0_rk)) delta_ij_plus_half = 0.0_rk + if (abs(delta_ij_plus_three_half) < epsilon(1.0_rk)) delta_ij_plus_three_half = 0.0_rk - r_L = delta_i_plus_half / (delta_i_minus_half + 1e-16_rk) - r_R = delta_i_plus_half / (delta_i_plus_three_half + 1e-16_rk) + r_L = delta_ij_plus_half / (delta_ij_minus_half + 1e-16_rk) + r_R = delta_ij_plus_half / (delta_ij_plus_three_half + 1e-16_rk) ! Inner term of alpha in Eq. 17 in Ref [1] #:if DIR == 'i' ! Xi variant, Eq 18a. Xi is used for the i-direction - tan_theta_i = abs(q(i, j + 1) - q(i, j - 1)) / (abs(q(i + 1, j) - q(i - 1, j)) + 1e-16_rk) - tan_theta_ip1 = abs(q(i + 1, j + 1) - q(i + 1, j - 1)) / (abs(q(i + 2, j) - q(i, j)) + 1e-16_rk) + tan_theta_ij = abs(q(i, j + 1) - q(i, j - 1)) / (abs(q(i + 1, j) - q(i - 1, j)) + 1e-16_rk) + tan_theta_ijp1 = abs(q(i + 1, j + 1) - q(i + 1, j - 1)) / (abs(q(i + 2, j) - q(i, j)) + 1e-16_rk) - ! phi_MN = [q(i,j+1), q(i+1,j+1), q(i,j), q(i+1,j), q(i,j-1), q(i+1,j+1)] + phi_MN = [q(i,j+1), q(i+1,j+1), q(i,j), q(i+1,j), q(i,j-1), q(i+1,j+1)] - min_phi_MN = huge(1.0_rk) - max_phi_MN = -huge(1.0_rk) - do n = j - 1, j + 1 - do m = i, i + 1 - if (q(m,n) < min_phi_MN) min_phi_MN = q(m,n) - if (q(m,n) > max_phi_MN) max_phi_MN = q(m,n) - end do - end do + ! min_phi_MN = huge(1.0_rk) + ! max_phi_MN = -huge(1.0_rk) + ! do n = j - 1, j + 1 + ! do m = i, i + 1 + ! if (q(m,n) < min_phi_MN) min_phi_MN = q(m,n) + ! if (q(m,n) > max_phi_MN) max_phi_MN = q(m,n) + ! end do + ! end do #:elif DIR == 'j' ! Eta variant, Eq 18b. Xi is used for the j-direction - tan_theta_i = abs(q(i + 1, j) - q(i - 1, j)) / (abs(q(i, j + 1) - q(i, j - 1)) + 1e-16_rk) - tan_theta_ip1 = abs(q(i + 1, j + 1) - q(i - 1, j + 1)) / (abs(q(i, j + 2) - q(i, j)) + 1e-16_rk) - - ! phi_MN = [q(i-1,j+1), q(i,j+1), q(i+1,j+1), q(i-1,j), q(i,j), q(i+1,j)] - - min_phi_MN = huge(1.0_rk) - max_phi_MN = -huge(1.0_rk) - do n = j, j + 1 - do m = i - 1, i + 1 - if (q(m,n) < min_phi_MN) min_phi_MN = q(m,n) - if (q(m,n) > max_phi_MN) max_phi_MN = q(m,n) - end do - end do + tan_theta_ij = abs(q(i + 1, j) - q(i - 1, j)) / (abs(q(i, j + 1) - q(i, j - 1)) + 1e-16_rk) + tan_theta_ijp1 = abs(q(i + 1, j + 1) - q(i - 1, j + 1)) / (abs(q(i, j + 2) - q(i, j)) + 1e-16_rk) + + phi_MN = [q(i-1,j+1), q(i,j+1), q(i+1,j+1), q(i-1,j), q(i,j), q(i+1,j)] + + ! min_phi_MN = huge(1.0_rk) + ! max_phi_MN = -huge(1.0_rk) + ! do n = j, j + 1 + ! do m = i - 1, i + 1 + ! if (q(m,n) < min_phi_MN) min_phi_MN = q(m,n) + ! if (q(m,n) > max_phi_MN) max_phi_MN = q(m,n) + ! end do + ! end do #:endif - if (delta_i_plus_half > 0.0_rk) then + min_phi_MN = minval(phi_MN) + max_phi_MN = maxval(phi_MN) + + if (delta_ij_plus_half > 0.0_rk) then phi_ML = max_phi_MN phi_MR = min_phi_MN else ! <= 0 @@ -473,33 +478,34 @@ subroutine interp_muscl_e_mlp_${DIR}$_edge_${ORDER}$_order(self, q, lbounds, edg end if ! Term inside the [..] in Eq 16a/17a in Ref [1] - alpha_L_term = (2.0_rk * max(1.0_rk, abs(r_L)) * abs((phi_ML - q(i, j)) / delta_i_plus_half)) / (1.0_rk + max(0.0_rk, tan_theta_i)) + alpha_L_term = (2.0_rk * max(1.0_rk, abs(r_L)) * abs((phi_ML - q(i, j)) / (delta_ij_plus_half + 1e-16_rk))) / (1.0_rk + max(0.0_rk, tan_theta_ij)) #:if DIR == 'i' ! Term inside the [..] in Eq 16b in Ref [1] - alpha_R_term = (2.0_rk * max(1.0_rk, abs(r_R)) * abs((q(i + 1, j) - phi_MR) / delta_i_plus_half)) / (1.0_rk + max(0.0_rk, tan_theta_ip1)) + alpha_R_term = (2.0_rk * max(1.0_rk, abs(r_R)) * abs((q(i + 1, j) - phi_MR) / (delta_ij_plus_half + 1e-16_rk))) / (1.0_rk + max(0.0_rk, tan_theta_ijp1)) #:elif DIR == 'j' ! Term inside the [..] in Eq 17b in Ref [1], there's a typo in the formula in the paper on the q(i + 1, j) term which should be q(i, j + 1) - alpha_R_term = (2.0_rk * max(1.0_rk, abs(r_R)) * abs((q(i, j + 1) - phi_MR) / delta_i_plus_half)) / (1.0_rk + max(0.0_rk, tan_theta_ip1)) + alpha_R_term = (2.0_rk * max(1.0_rk, abs(r_R)) * abs((q(i, j + 1) - phi_MR) / (delta_ij_plus_half + 1e-16_rk))) / (1.0_rk + max(0.0_rk, tan_theta_ijp1)) #:endif ! This is the g(x) = max(1, min(2, alpha)) function alpha_L = max(1.0_rk, min(2.0_rk, alpha_L_term)) alpha_R = max(1.0_rk, min(2.0_rk, alpha_R_term)) + ! write(*, '(a, 10(es16.6))') '${DIR}$, ', alpha_L_term, alpha_R_term, alpha_L, alpha_R ! Phi L, e.g. left state - edge_values(1, i, j) = q(i, j) + 0.5_rk * max(0.0_rk, min(alpha_L * delta_i_plus_half, & - alpha_L * delta_i_minus_half, & - beta_L_delta_i_minus_half)) + edge_values(1, i, j) = q(i, j) + 0.5_rk * max(0.0_rk, min(alpha_L * delta_ij_plus_half, & + alpha_L * delta_ij_minus_half, & + beta_L_delta_ij_minus_half)) ! Phi R, e.g. right state #:if DIR == 'i' - edge_values(2, i, j) = q(i + 1, j) - 0.5_rk * max(0.0_rk, min(alpha_R * delta_i_plus_half, & - alpha_R * delta_i_plus_three_half, & - beta_R_delta_i_plus_three_half)) + edge_values(2, i, j) = q(i + 1, j) - 0.5_rk * max(0.0_rk, min(alpha_R * delta_ij_plus_half, & + alpha_R * delta_ij_plus_three_half, & + beta_R_delta_ij_plus_three_half)) #:elif DIR == 'j' - edge_values(2, i, j) = q(i, j + 1) - 0.5_rk * max(0.0_rk, min(alpha_R * delta_i_plus_half, & - alpha_R * delta_i_plus_three_half, & - beta_R_delta_i_plus_three_half)) + edge_values(2, i, j) = q(i, j + 1) - 0.5_rk * max(0.0_rk, min(alpha_R * delta_ij_plus_half, & + alpha_R * delta_ij_plus_three_half, & + beta_R_delta_ij_plus_three_half)) #:endif end select diff --git a/src/lib/spatial_reconstruction/muscl_mlp.fypp b/src/lib/spatial_reconstruction/muscl_mlp.fypp index 5fa5469..e541da9 100644 --- a/src/lib/spatial_reconstruction/muscl_mlp.fypp +++ b/src/lib/spatial_reconstruction/muscl_mlp.fypp @@ -25,8 +25,10 @@ #ifdef __SIMD_ALIGN_OMP__ #define __INTERP_ALIGN__ aligned(q, edge_values:__ALIGNBYTES__) +#define __CONT_ALIGN__ aligned(rho, u, v, p, d_bar_rho, d_bar_u, d_bar_v, d_bar_p:__ALIGNBYTES__) #else #define __INTERP_ALIGN__ +#define __CONT_ALIGN__ #endif module mod_muscl_mlp @@ -47,7 +49,6 @@ module mod_muscl_mlp public :: muscl_mlp_t, new_muscl_mlp type, extends(muscl_interpolation_t) :: muscl_mlp_t - integer(ik), dimension(:,:), allocatable :: continuity_sensor !< (i,j); flag for continuous and linear/non-linear discontinuity regions contains procedure, public :: initialize => init_muscl_mlp procedure, public :: interpolate_edge_values @@ -64,20 +65,35 @@ module mod_muscl_mlp contains -function new_muscl_mlp(limiter) result(interpolator) +function new_muscl_mlp(limiter, order) result(interpolator) type(muscl_mlp_t), pointer :: interpolator character(len=*), intent(in) :: limiter + integer(ik), intent(in) :: order allocate(interpolator) - interpolator%limiter_name = trim(limiter) - interpolator%order = 2 + write(interpolator%limiter_name, '(a, i0)') 'MLP', order + interpolator%order = order + end function subroutine init_muscl_mlp(self, limiter) class(muscl_mlp_t), intent(inout) :: self character(len=*), intent(in) :: limiter - self%limiter_name = trim(limiter) - self%order = 2 + + + select case(trim(limiter)) + case('MLP3') + self%limiter_name = 'MLP3' + self%order = 3 + case('MLP5') + self%limiter_name = 'MLP5' + self%order = 5 + case default + call error_msg(module='mod_muscl_mlp', class='muscl_mlp_t', procedure='init_muscl_mlp', & + message="Unknown limiter type: '" // trim(limiter) // "'", & + file_name=__FILE__, line_number=__LINE__) + end select + end subroutine init_muscl_mlp subroutine interpolate_edge_values(self, q, lbounds, i_edges, j_edges) @@ -90,29 +106,37 @@ subroutine interpolate_edge_values(self, q, lbounds, i_edges, j_edges) real(rk), dimension(:, :, :), allocatable, intent(out) :: i_edges real(rk), dimension(:, :, :), allocatable, intent(out) :: j_edges - ! #:for DIR in DIRECTIONS - ! call self%interp_muscl_mlp_${DIR}$_edge(q, lbounds, ${DIR}$_edges) - ! #:endfor + select case(self%order) + #:for ORDER in SPATIAL_ORDER + case(${ORDER}$) + #:for DIR in DIRECTIONS + call self%interp_muscl_mlp_${DIR}$_edge_${ORDER}$_order(q, lbounds, ${DIR}$_edges) + #:endfor + #:endfor + case default + call error_msg(module='mod_muscl_mlp', class='muscl_mlp_t', procedure='interpolate_edge_values', & + message="Unknown limiter type: '" // trim(self%limiter_name) // "'", & + file_name=__FILE__, line_number=__LINE__) + end select end subroutine interpolate_edge_values subroutine distinguish_continuous_regions(self, rho, u, v, p, lbounds) - !< For plain-jane TVD2, there is no need to distinguish the regions which are continuous or linear/non-linear discontinuity regions. This - !< is required by the parent class, but is really only implemented in the e-MLP classes. + !< Find the continuous and discontinuous regions class(muscl_mlp_t), intent(inout) :: self integer(ik), dimension(2), intent(in) :: lbounds real(rk), dimension(lbounds(1):, lbounds(2):), contiguous, intent(in) :: rho !< (i,j); density - real(rk), dimension(lbounds(1):, lbounds(2):), contiguous, intent(in) :: u !< (i,j); density - real(rk), dimension(lbounds(1):, lbounds(2):), contiguous, intent(in) :: v !< (i,j); density - real(rk), dimension(lbounds(1):, lbounds(2):), contiguous, intent(in) :: p !< (i,j); density + real(rk), dimension(lbounds(1):, lbounds(2):), contiguous, intent(in) :: u !< (i,j); x-velocity + real(rk), dimension(lbounds(1):, lbounds(2):), contiguous, intent(in) :: v !< (i,j); y-velocity + real(rk), dimension(lbounds(1):, lbounds(2):), contiguous, intent(in) :: p !< (i,j); pressure end subroutine distinguish_continuous_regions subroutine finalize(self) !< Finalizer routine to clean up allocated data type(muscl_mlp_t), intent(inout) :: self - if (allocated(self%continuity_sensor)) deallocate(self%continuity_sensor) end subroutine finalize + #:for DIR in DIRECTIONS #:for ORDER in SPATIAL_ORDER subroutine interp_muscl_mlp_${DIR}$_edge_${ORDER}$_order(self, q, lbounds, edge_values) @@ -125,17 +149,27 @@ subroutine interp_muscl_mlp_${DIR}$_edge_${ORDER}$_order(self, q, lbounds, edge_ real(rk), dimension(:, :, :), allocatable, intent(out) :: edge_values !<((L,R), i, j); L/R state for each edge - integer(ik) :: i, j + integer(ik) :: i, j, m, n integer(ik) :: ilo, ihi, jlo, jhi - real(rk) :: delta_i_minus_half, delta_i_plus_half, delta_i_plus_three_half - real(rk) :: r_R, r_L, phi_limit + real(rk) :: delta_ij_minus_half !< Delta Phi_{i-1/2, j} = q(i, j) - q(i - 1, j) (or j-1/2 depending on direction) + real(rk) :: delta_ij_plus_half !< Delta Phi_{i+1/2, j} = q(i + 1, j) - q(i, j) (or j+1/2 depending on direction) + real(rk) :: delta_ij_plus_three_half !< Delta Phi_{i+3/2, j} = q(i + 2, j) - q(i + 1, j) (or j+3/2 depending on direction) + + + real(rk) :: r_L_ij, r_R_ij_p1, r_L_ij_p1, r_L_ij_m1, r_R_ij, r_R_ij_p2 + + real(rk) :: tan_theta_ij !< the tan(theta), i term + real(rk) :: tan_theta_ij_p1 !< the tan(theta), i+1 term + real(rk) :: beta_L, beta_R + real(rk) :: alpha_L, alpha_R + real(rk) :: alpha_L_term, alpha_R_term !dir$ assume_aligned q: __ALIGNBYTES__ !dir$ assume_aligned edge_values: __ALIGNBYTES__ - if (enable_debug_print) call debug_print('Running muscl_tvd2_t%interp_muscl_mlp_${DIR}$_edge()', __FILE__, __LINE__) + if (enable_debug_print) call debug_print('Running muscl_tvd2_t%interp_muscl_mlp_${DIR}$_edge_${ORDER}$_order()', __FILE__, __LINE__) ! Because this is an edge-based array, the indexing is # cells + 1, and ! only goes into the single nearest ghost layer @@ -145,9 +179,9 @@ subroutine interp_muscl_mlp_${DIR}$_edge_${ORDER}$_order(self, q, lbounds, edge_ jhi = ubound(q, dim=2) - n_ghost_layers #:if DIR == 'i' - allocate(edge_values(1:2, ilo-1:ihi, jlo:jhi)) + allocate(edge_values(1:2, ilo-1:ihi, jlo:jhi)) #:elif DIR == 'j' - allocate(edge_values(1:2, ilo:ihi, jlo-1:jhi)) + allocate(edge_values(1:2, ilo:ihi, jlo-1:jhi)) #:endif edge_values = 0.0_rk @@ -164,78 +198,118 @@ subroutine interp_muscl_mlp_${DIR}$_edge_${ORDER}$_order(self, q, lbounds, edge_ ! ! This is the numbering convention that this module uses - !$omp parallel default(none), & - !$omp firstprivate(ilo, ihi, jlo, jhi) & - !$omp private(i, j, phi_limit, delta_i_minus_half, delta_i_plus_half, delta_i_plus_three_half) & - !$omp private(r_L, r_R) & - !$omp shared(q, edge_values) - !$omp do - do j = jlo, jhi - !$omp simd __INTERP_ALIGN__ - !dir$ vector aligned - do i = ilo, ihi - - #:if DIR == 'i' - delta_i_minus_half = q(i, j) - q(i - 1, j) - delta_i_plus_half = q(i + 1, j) - q(i, j) - delta_i_plus_three_half = q(i + 2, j) - q(i + 1, j) - #:elif DIR == 'j' - delta_i_minus_half = q(i, j) - q(i, j - 1) - delta_i_plus_half = q(i, j + 1) - q(i, j) - delta_i_plus_three_half = q(i, j + 2) - q(i, j + 1) - #:endif - - if (abs(delta_i_minus_half) < epsilon(1.0_rk)) delta_i_minus_half = 0.0_rk - if (abs(delta_i_plus_half) < epsilon(1.0_rk)) delta_i_plus_half = 0.0_rk - if (abs(delta_i_plus_three_half) < epsilon(1.0_rk)) delta_i_plus_three_half = 0.0_rk - - r_L = delta_i_plus_half / (delta_i_minus_half + 1e-16_rk) - r_R = delta_i_plus_half / (delta_i_plus_three_half + 1e-16_rk) - - - ! #:if ORDER == 3 - ! ! 3rd order interpolation - ! beta_L = (1.0_rk + 2.0_rk * r_L_i) / 3.0_rk - ! beta_R = (1.0_rk + 2.0_rk * r_R_ip1) / 3.0_rk - ! #:elif ORDER == 5 - ! ! 5th order interpolation - ! beta_L = ((-2.0_rk / r_L_im1) + 11.0_rk + (24.0_rk * r_L_i) - (3.0_rk * r_L_i * r_L_ip1)) / 30.0_rk - ! beta_R = ((-2.0_rk / r_R_ip2) + 11.0_rk + (24.0_rk * r_R_ip1) - (3.0_rk * r_R_ip1 * r_R_i)) / 30.0_rk - ! #:endif - - - ! ! Apply the limiter - ! #:if LIMITER == 'minmod' - ! phi_limit = max(0.0_rk, min(r_L, 1.0_rk)) - ! #:elif LIMITER == 'superbee' - ! phi_limit = max(0.0_rk, min(2.0_rk * r_L, 1.0_rk), min(r_L, 2.0_rk)) - ! #:elif LIMITER == 'vanleer' - ! phi_limit = (r_L + abs(r_L)) / (1.0_rk + abs(r_L)) - ! #:endif - - ! ! Phi L, e.g. left state - ! edge_values(1, i, j) = q(i, j) + 0.5_rk * phi_limit * delta_i_minus_half - - ! ! Apply the limiter - ! #:if LIMITER == 'minmod' - ! phi_limit = max(0.0_rk, min(r_R, 1.0_rk)) - ! #:elif LIMITER == 'superbee' - ! phi_limit = max(0.0_rk, min(2.0_rk * r_R, 1.0_rk), min(r_R, 2.0_rk)) - ! #:elif LIMITER == 'vanleer' - ! phi_limit = (r_R + abs(r_R)) / (1.0_rk + abs(r_R)) - ! #:endif - - ! ! Phi R, e.g. right state - ! #:if DIR == 'i' - ! edge_values(2, i, j) = q(i + 1, j) - 0.5_rk * phi_limit * delta_i_plus_three_half - ! #:elif DIR == 'j' - ! edge_values(2, i, j) = q(i, j + 1) - 0.5_rk * phi_limit * delta_i_plus_three_half - ! #:endif - - end do - end do - !$omp end do - !$omp end parallel + ! !$omp parallel default(none), & + ! !$omp firstprivate(ilo, ihi, jlo, jhi) & + ! !$omp private(i, j) & + ! !$omp private(delta_ij_minus_half, delta_ij_plus_half, delta_ij_plus_three_half ) & + ! !$omp private(r_L_ij, r_R_ij_p1, r_L_ij_p1, r_L_ij_m1, r_R_ij, r_R_ij_p2) & + ! !$omp private(tan_theta_ij, tan_theta_ij_p1, beta_L, beta_R, alpha_L, alpha_R, alpha_L_term, alpha_R_term) & + ! !$omp shared(q, edge_values, self) + ! !$omp do + ! do j = jlo, jhi + ! !! $omp simd __INTERP_ALIGN__ + ! !! dir$ vector aligned + ! do i = ilo, ihi + + ! #:if DIR == 'i' + ! delta_ij_minus_half = q(i, j) - q(i - 1, j) + ! delta_ij_plus_half = q(i + 1, j) - q(i, j) + ! delta_ij_plus_three_half = q(i + 2, j) - q(i + 1, j) + ! #:elif DIR == 'j' + ! delta_ij_minus_half = q(i, j) - q(i, j - 1) + ! delta_ij_plus_half = q(i, j + 1) - q(i, j) + ! delta_ij_plus_three_half = q(i, j + 2) - q(i, j + 1) + ! #:endif + + ! if (abs(delta_ij_minus_half) < epsilon(1.0_rk)) delta_ij_minus_half = 0.0_rk + ! if (abs(delta_ij_plus_half) < epsilon(1.0_rk)) delta_ij_plus_half = 0.0_rk + ! if (abs(delta_ij_plus_three_half) < epsilon(1.0_rk)) delta_ij_plus_three_half = 0.0_rk + + ! r_L_ij = delta_ij_plus_half / (delta_ij_minus_half + 1e-16_rk) + ! r_R_ij_p1 = delta_ij_plus_half / (delta_ij_plus_three_half + 1e-16_rk) + + ! #:if ORDER == 3 + ! ! 3rd order interpolation + ! ! Eq 65a in Ref [1] + ! beta_L = (1.0_rk + 2.0_rk * r_L_ij) / 3.0_rk + + ! ! Eq 65b in Ref [1] + ! beta_R = (1.0_rk + 2.0_rk * r_R_ij_p1) / 3.0_rk + ! #:elif ORDER == 5 + ! error stop "5th order not working yet" + ! r_L_ij_p1 = 0 + ! r_L_ij_m1 = 0 + ! r_R_ij = 0 + ! r_R_ij_p2 = 0 + + ! ! Eq 23a in Ref [1] + ! beta_L = ((-2.0_rk / r_L_ij_m1) & + ! + 11.0_rk & + ! + (24.0_rk * r_L_ij) & + ! - (3.0_rk * r_L_ij * r_L_ij_p1)) / 30.0_rk + + ! ! Eq 23b in Ref [1] + ! beta_R = ((-2.0_rk / r_R_ij_p2) & + ! + 11.0_rk & + ! + (24.0_rk * r_R_ij_p1) & + ! - (3.0_rk * r_R_ij_p1 * r_R_ij)) / 30.0_rk + ! #:endif + + ! ! Inner term of alpha in Eq. 17 in Ref [1] + ! #:if DIR == 'i' + ! ! Xi variant, Eq 18a. Xi is used for the i-direction + ! tan_theta_ij = abs(q(i, j + 1) - q(i, j - 1)) / (abs(q(i + 1, j) - q(i - 1, j)) + 1e-16_rk) + ! tan_theta_ij_p1 = abs(q(i + 1, j + 1) - q(i + 1, j - 1)) / (abs(q(i + 2, j) - q(i, j)) + 1e-16_rk) + ! #:elif DIR == 'j' + ! ! Eta variant, Eq 18b. Xi is used for the j-direction + ! tan_theta_ij = abs(q(i + 1, j) - q(i - 1, j)) / (abs(q(i, j + 1) - q(i, j - 1)) + 1e-16_rk) + ! tan_theta_ij_p1 = abs(q(i + 1, j + 1) - q(i - 1, j + 1)) / (abs(q(i, j + 2) - q(i, j)) + 1e-16_rk) + ! #:endif + + ! #:if DIR == 'i' + ! ! Term inside the [..] in Eq 64a + ! alpha_L_term = (2.0_rk * max(1.0_rk, r_L_ij) * & + ! (1.0_rk + max(0.0_rk, (tan_theta_ij_p1 / r_R_ij_p1)))) / & + ! (1.0_rk + tan_theta_ij) + + ! alpha_R_term = (2.0_rk * max(1.0_rk, r_R_ij_p1) * & + ! (1.0_rk + max(0.0_rk, (tan_theta_ij / r_L_ij)))) / & + ! (1.0_rk + tan_theta_ij_p1) + ! #:elif DIR == 'j' + ! ! Term inside the [..] in Eq 64b + ! alpha_L_term = (2.0_rk * max(1.0_rk, r_L_ij) * & + ! (1.0_rk + max(0.0_rk, (tan_theta_ij_p1 / r_R_ij_p1)))) / & + ! (1.0_rk + tan_theta_ij) + + ! alpha_R_term = (2.0_rk * max(1.0_rk, r_R_ij_p1) * & + ! (1.0_rk + max(0.0_rk, (tan_theta_ij / r_L_ij)))) / & + ! (1.0_rk + tan_theta_ij_p1) + ! #:endif + + ! ! This is the g(x) = max(1, min(2, alpha)) function + ! alpha_L = max(1.0_rk, min(2.0_rk, alpha_L_term)) + ! alpha_R = max(1.0_rk, min(2.0_rk, alpha_R_term)) + ! ! write(*, '(a, 10(es16.6))') '${DIR}$, ', alpha_L_term, alpha_R_term, alpha_L, alpha_R + + ! ! Phi L, e.g. left state + ! edge_values(1, i, j) = q(i, j) + 0.5_rk * max(0.0_rk, min(alpha_L * delta_ij_plus_half, & + ! alpha_L * delta_ij_minus_half, & + ! beta_L_delta_ij_minus_half)) + ! ! Phi R, e.g. right state + ! #:if DIR == 'i' + ! edge_values(2, i, j) = q(i + 1, j) - 0.5_rk * max(0.0_rk, min(alpha_R * r_R_ip1, & + ! alpha_R, & + ! beta_R)) * delta_ij_plus_three_half + ! #:elif DIR == 'j' + ! edge_values(2, i, j) = q(i, j + 1) - 0.5_rk * max(0.0_rk, min(alpha_R * r_R_ip1, & + ! alpha_R, & + ! beta_R)) * delta_ij_plus_three_half + ! #:endif + + ! end do + ! end do + ! !$omp end do + ! !$omp end parallel end subroutine interp_muscl_mlp_${DIR}$_edge_${ORDER}$_order