diff --git a/docs/documentation/thermochemistry.md b/docs/documentation/thermochemistry.md index 7398d0593..12a3bac0f 100644 --- a/docs/documentation/thermochemistry.md +++ b/docs/documentation/thermochemistry.md @@ -31,7 +31,16 @@ to parse a mechanism does not imply that the generator supports every feature in The generated module provides species metadata, caloric and ideal-gas properties, temperature inversion, net production rates, fused creation/destruction rates, -mixture viscosity and thermal conductivity, and species diffusivities. MFC still +mixture viscosity and thermal conductivity, and species diffusivities. Routines +shaped for MFC's callers return what a call site needs from one state: +`get_mixture_caloric_state` (species cp/R and mixture cp, cv and energy from one +pass over the NASA7 fits), `get_mixavg_transport_state` (molecular weight, mole +fractions, mixture-averaged diffusivities and conductivity, sharing the composition +work) and `get_species_enthalpies_mass`. They live in the companion module +`m_thermochem_state`, generated into the same file, so `m_thermochem`, which nearly +every MFC file reads, is unchanged; nvfortran 24.7 and older crash compiling unrelated +files when it grows by these routines. They reuse the separate routines' +arithmetic; optimized builds may still differ from separate calls at roundoff. MFC still owns reaction time integration, including alpha-QSS, and spatial transport discretization. The generator emits single- or double-precision routines and the selected CPU, OpenACC or OpenMP annotations. Nonchemistry builds retain the existing diff --git a/src/common/m_chemistry.fpp b/src/common/m_chemistry.fpp index 687b2a118..2735c712a 100644 --- a/src/common/m_chemistry.fpp +++ b/src/common/m_chemistry.fpp @@ -14,7 +14,7 @@ module m_chemistry & get_species_mass_diffusivities_mixavg, gas_constant, get_mixture_molecular_weight, get_mixture_energy_mass, & & get_mixture_thermal_conductivity_mixavg, get_species_enthalpies_rt, get_mixture_viscosity_mixavg, & & get_mixture_specific_heat_cp_mass, get_mixture_enthalpy_mass - + use m_thermochem_state, only: get_mixavg_transport_state, get_species_enthalpies_mass use m_global_parameters implicit none @@ -385,14 +385,6 @@ contains & - eqn_idx%species%beg + 1)) end do - ! Calculate molecular weights and mole fractions - call get_mixture_molecular_weight(Ys_L, MW_L) - call get_mixture_molecular_weight(Ys_R, MW_R) - MW_cell = 0.5_wp*(MW_L + MW_R) - - call get_mole_fractions(MW_L, Ys_L, Xs_L) - call get_mole_fractions(MW_R, Ys_R, Xs_R) - P_L = q_prim_qp(eqn_idx%E)%sf(x, y, z) P_R = q_prim_qp(eqn_idx%E)%sf(x + offsets(1), y + offsets(2), z + offsets(3)) @@ -405,23 +397,17 @@ contains rho_cell = 0.5_wp*(rho_L + rho_R) dT_dxi = (T_R - T_L)/grid_spacing - ! Get transport properties - call get_species_mass_diffusivities_mixavg(P_L, T_L, Ys_L, mass_diffusivities_mixavg1) - call get_species_mass_diffusivities_mixavg(P_R, T_R, Ys_R, mass_diffusivities_mixavg2) - - call get_mixture_thermal_conductivity_mixavg(T_L, Ys_L, lambda_L) - call get_mixture_thermal_conductivity_mixavg(T_R, Ys_R, lambda_R) + ! Transport state and species enthalpies [J/kg] of each side + call get_mixavg_transport_state(P_L, T_L, Ys_L, MW_L, Xs_L, mass_diffusivities_mixavg1, lambda_L) + call get_mixavg_transport_state(P_R, T_R, Ys_R, MW_R, Xs_R, mass_diffusivities_mixavg2, lambda_R) + MW_cell = 0.5_wp*(MW_L + MW_R) - call get_species_enthalpies_rt(T_L, h_l) - call get_species_enthalpies_rt(T_R, h_r) + call get_species_enthalpies_mass(T_L, h_l) + call get_species_enthalpies_mass(T_R, h_r) ! Calculate species properties and gradients $:GPU_LOOP(parallelism='[seq]') do i = eqn_idx%species%beg, eqn_idx%species%end - h_l(i - eqn_idx%species%beg + 1) = h_l(i - eqn_idx%species%beg + 1) & - & *gas_constant*T_L/molecular_weights(i - eqn_idx%species%beg + 1) - h_r(i - eqn_idx%species%beg + 1) = h_r(i - eqn_idx%species%beg + 1) & - & *gas_constant*T_R/molecular_weights(i - eqn_idx%species%beg + 1) Xs_cell(i - eqn_idx%species%beg + 1) = 0.5_wp*(Xs_L(i - eqn_idx%species%beg + 1) + Xs_R(i & & - eqn_idx%species%beg + 1)) h_k(i - eqn_idx%species%beg + 1) = 0.5_wp*(h_l(i - eqn_idx%species%beg + 1) + h_r(i & diff --git a/src/simulation/m_cbc.fpp b/src/simulation/m_cbc.fpp index 7bddbdc31..705cca23c 100644 --- a/src/simulation/m_cbc.fpp +++ b/src/simulation/m_cbc.fpp @@ -15,9 +15,9 @@ module m_cbc use m_compute_cbc use m_boundary_primitives, only: f_vel_ramp use m_constants, only: riemann_solver_hll, model_eqns_gamma_law, recon_type_weno, recon_type_muscl - use m_thermochem, only: get_mixture_energy_mass, get_mixture_specific_heat_cv_mass, get_mixture_specific_heat_cp_mass, & - & gas_constant, get_mixture_molecular_weight, get_species_enthalpies_rt, molecular_weights, get_species_specific_heats_r, & + use m_thermochem, only: gas_constant, get_mixture_molecular_weight, get_species_enthalpies_rt, molecular_weights, & & get_mole_fractions + use m_thermochem_state, only: get_mixture_caloric_state implicit none @@ -646,18 +646,15 @@ contains call get_mixture_molecular_weight(Ys, Mw) R_gas = gas_constant/Mw T = pres/rho/R_gas - call get_mixture_specific_heat_cp_mass(T, Ys, Cp) - call get_mixture_energy_mass(T, Ys, e_mix) + call get_mixture_caloric_state(T, Ys, Cp_i, Cp, Cv, e_mix) E = rho*e_mix + 5.e-1_wp*rho*vel_K_sum if (chem_params%gamma_method == 1) then !> gamma_method = 1: Ref. Section 2.3.1 Formulation of doi:10.7907/ZKW8-ES97. call get_mole_fractions(Mw, Ys, Xs) - call get_species_specific_heats_r(T, Cp_i) Gamma_i(1:num_species) = Cp_i(1:num_species)/(Cp_i(1:num_species) - 1.0_wp) gamma = sum(Xs(1:num_species)/(Gamma_i(1:num_species) - 1.0_wp)) else if (chem_params%gamma_method == 2) then !> gamma_method = 2: c_p / c_v where c_p, c_v are specific heats. - call get_mixture_specific_heat_cv_mass(T, Ys, Cv) gamma = 1.0_wp/(Cp/Cv - 1.0_wp) end if end if diff --git a/src/simulation/m_riemann_solver_hll.fpp b/src/simulation/m_riemann_solver_hll.fpp index 969d6996a..b9ea17387 100644 --- a/src/simulation/m_riemann_solver_hll.fpp +++ b/src/simulation/m_riemann_solver_hll.fpp @@ -15,9 +15,8 @@ module m_riemann_solver_hll use m_constants, only: riemann_solver_hll, riemann_solver_hllc, riemann_solver_lax_friedrichs, avg_state_roe, & & avg_state_arithmetic, wave_speeds_direct, wave_speeds_pressure use m_chemistry - use m_thermochem, only: gas_constant, get_mixture_molecular_weight, get_mixture_specific_heat_cv_mass, & - & get_mixture_energy_mass, get_species_specific_heats_r, get_species_enthalpies_rt, get_mixture_specific_heat_cp_mass, & - & molecular_weights + use m_thermochem, only: gas_constant, get_mixture_molecular_weight, get_species_enthalpies_rt, molecular_weights + use m_thermochem_state, only: get_mixture_caloric_state use m_riemann_state implicit none @@ -224,8 +223,8 @@ contains T_L = pres_L/rho_L/R_gas_L T_R = pres_R/rho_R/R_gas_R - call get_species_specific_heats_r(T_L, Cp_iL) - call get_species_specific_heats_r(T_R, Cp_iR) + call get_mixture_caloric_state(T_L, Ys_L, Cp_iL, Cp_L, Cv_L, E_L) + call get_mixture_caloric_state(T_R, Ys_R, Cp_iR, Cp_R, Cv_R, E_R) if (chem_params%gamma_method == 1) then ! gamma_method = 1: Ref. Section 2.3.1 Formulation of doi:10.7907/ZKW8-ES97. @@ -236,20 +235,12 @@ contains gamma_R = sum(Xs_R(1:num_species)/(Gamma_iR(1:num_species) - 1.0_wp)) else if (chem_params%gamma_method == 2) then ! gamma_method = 2: c_p / c_v where c_p, c_v are specific heats. - call get_mixture_specific_heat_cp_mass(T_L, Ys_L, Cp_L) - call get_mixture_specific_heat_cp_mass(T_R, Ys_R, Cp_R) - call get_mixture_specific_heat_cv_mass(T_L, Ys_L, Cv_L) - call get_mixture_specific_heat_cv_mass(T_R, Ys_R, Cv_R) - Gamm_L = Cp_L/Cv_L gamma_L = 1.0_wp/(Gamm_L - 1.0_wp) Gamm_R = Cp_R/Cv_R gamma_R = 1.0_wp/(Gamm_R - 1.0_wp) end if - call get_mixture_energy_mass(T_L, Ys_L, E_L) - call get_mixture_energy_mass(T_R, Ys_R, E_R) - E_L = rho_L*E_L + 5.e-1*rho_L*vel_L_rms E_R = rho_R*E_R + 5.e-1*rho_R*vel_R_rms H_L = (E_L + pres_L)/rho_L diff --git a/src/simulation/m_riemann_solver_hllc.fpp b/src/simulation/m_riemann_solver_hllc.fpp index d96c8da7e..6d49f2931 100644 --- a/src/simulation/m_riemann_solver_hllc.fpp +++ b/src/simulation/m_riemann_solver_hllc.fpp @@ -18,9 +18,8 @@ module m_riemann_solver_hllc use m_bubbles_EE use m_surface_tension use m_chemistry - use m_thermochem, only: gas_constant, get_mixture_molecular_weight, get_mixture_specific_heat_cv_mass, & - & get_mixture_energy_mass, get_species_specific_heats_r, get_species_enthalpies_rt, get_mixture_specific_heat_cp_mass, & - & molecular_weights + use m_thermochem, only: gas_constant, get_mixture_molecular_weight, get_species_enthalpies_rt, molecular_weights + use m_thermochem_state, only: get_mixture_caloric_state use m_riemann_state implicit none @@ -1020,8 +1019,8 @@ contains T_L = pres_L/rho_L/R_gas_L T_R = pres_R/rho_R/R_gas_R - call get_species_specific_heats_r(T_L, Cp_iL) - call get_species_specific_heats_r(T_R, Cp_iR) + call get_mixture_caloric_state(T_L, Ys_L, Cp_iL, Cp_L, Cv_L, E_L) + call get_mixture_caloric_state(T_R, Ys_R, Cp_iR, Cp_R, Cv_R, E_R) if (chem_params%gamma_method == 1) then !> gamma_method = 1: Ref. Section 2.3.1 Formulation of doi:10.7907/ZKW8-ES97. @@ -1032,18 +1031,10 @@ contains gamma_R = sum(Xs_R(1:num_species)/(Gamma_iR(1:num_species) - 1.0_wp)) else if (chem_params%gamma_method == 2) then !> gamma_method = 2: c_p / c_v where c_p, c_v are specific heats. - call get_mixture_specific_heat_cp_mass(T_L, Ys_L, Cp_L) - call get_mixture_specific_heat_cp_mass(T_R, Ys_R, Cp_R) - call get_mixture_specific_heat_cv_mass(T_L, Ys_L, Cv_L) - call get_mixture_specific_heat_cv_mass(T_R, Ys_R, Cv_R) - Gamm_L = Cp_L/Cv_L; Gamm_R = Cp_R/Cv_R gamma_L = 1.0_wp/(Gamm_L - 1.0_wp); gamma_R = 1.0_wp/(Gamm_R - 1.0_wp) end if - call get_mixture_energy_mass(T_L, Ys_L, E_L) - call get_mixture_energy_mass(T_R, Ys_R, E_R) - E_L = rho_L*E_L + 5.e-1*rho_L*vel_L_rms E_R = rho_R*E_R + 5.e-1*rho_R*vel_R_rms H_L = (E_L + pres_L)/rho_L diff --git a/src/simulation/m_riemann_solver_lf.fpp b/src/simulation/m_riemann_solver_lf.fpp index 82dd93059..5672d955a 100644 --- a/src/simulation/m_riemann_solver_lf.fpp +++ b/src/simulation/m_riemann_solver_lf.fpp @@ -12,8 +12,8 @@ module m_riemann_solver_lf use m_global_parameters use m_variables_conversion use m_constants, only: riemann_solver_hll, riemann_solver_hllc, riemann_solver_lax_friedrichs - use m_thermochem, only: gas_constant, get_mixture_molecular_weight, get_mixture_specific_heat_cv_mass, & - & get_mixture_energy_mass, get_species_specific_heats_r, get_mixture_specific_heat_cp_mass, molecular_weights + use m_thermochem, only: gas_constant, get_mixture_molecular_weight, molecular_weights + use m_thermochem_state, only: get_mixture_caloric_state use m_riemann_state implicit none @@ -178,8 +178,8 @@ contains T_L = pres_L/rho_L/R_gas_L T_R = pres_R/rho_R/R_gas_R - call get_species_specific_heats_r(T_L, Cp_iL) - call get_species_specific_heats_r(T_R, Cp_iR) + call get_mixture_caloric_state(T_L, Ys_L, Cp_iL, Cp_L, Cv_L, E_L) + call get_mixture_caloric_state(T_R, Ys_R, Cp_iR, Cp_R, Cv_R, E_R) if (chem_params%gamma_method == 1) then ! gamma_method = 1: Ref. Section 2.3.1 Formulation of doi:10.7907/ZKW8-ES97. @@ -190,20 +190,12 @@ contains gamma_R = sum(Xs_R(1:num_species)/(Gamma_iR(1:num_species) - 1.0_wp)) else if (chem_params%gamma_method == 2) then ! gamma_method = 2: c_p / c_v where c_p, c_v are specific heats. - call get_mixture_specific_heat_cp_mass(T_L, Ys_L, Cp_L) - call get_mixture_specific_heat_cp_mass(T_R, Ys_R, Cp_R) - call get_mixture_specific_heat_cv_mass(T_L, Ys_L, Cv_L) - call get_mixture_specific_heat_cv_mass(T_R, Ys_R, Cv_R) - Gamm_L = Cp_L/Cv_L gamma_L = 1.0_wp/(Gamm_L - 1.0_wp) Gamm_R = Cp_R/Cv_R gamma_R = 1.0_wp/(Gamm_R - 1.0_wp) end if - call get_mixture_energy_mass(T_L, Ys_L, E_L) - call get_mixture_energy_mass(T_R, Ys_R, E_R) - E_L = rho_L*E_L + 5.e-1*rho_L*vel_L_rms E_R = rho_R*E_R + 5.e-1*rho_R*vel_R_rms else diff --git a/toolchain/mfc/test_thermochem.py b/toolchain/mfc/test_thermochem.py index 17a0f5f34..d31b50960 100644 --- a/toolchain/mfc/test_thermochem.py +++ b/toolchain/mfc/test_thermochem.py @@ -223,3 +223,56 @@ def test_zero_concentration_falloff(tmp_path, precision): for i, reaction in enumerate(gas.reactions()): if isinstance(reaction.rate, ct.FalloffRate): assert rates[i] == 0 + + +FUSED_DRIVER = """ +program fused + use m_thermochem + use m_thermochem_state + implicit none + integer :: ierr + real(KIND) :: t, pressure, cp, cv, e, lambda, mw, cp_f, cv_f, e_f, lambda_f, mw_f + real(KIND), dimension(num_species) :: y, cp0_r, cp0_r_f, x, x_f, diffusion, diffusion_f, h_rt, h + do + read(*,*,iostat=ierr) t, pressure, y + if (ierr /= 0) exit + call get_species_specific_heats_r(t, cp0_r) + call get_mixture_specific_heat_cp_mass(t, y, cp) + call get_mixture_specific_heat_cv_mass(t, y, cv) + call get_mixture_energy_mass(t, y, e) + call get_mixture_caloric_state(t, y, cp0_r_f, cp_f, cv_f, e_f) + call get_mixture_molecular_weight(y, mw) + call get_mole_fractions(mw, y, x) + call get_species_mass_diffusivities_mixavg(pressure, t, y, diffusion) + call get_mixture_thermal_conductivity_mixavg(t, y, lambda) + call get_mixavg_transport_state(pressure, t, y, mw_f, x_f, diffusion_f, lambda_f) + call get_species_enthalpies_rt(t, h_rt) + call get_species_enthalpies_mass(t, h) + ! Fused routines share the separate routines' arithmetic, so they agree exactly. + if (any(cp0_r /= cp0_r_f) .or. cp /= cp_f .or. cv /= cv_f .or. e /= e_f) stop 1 + if (mw /= mw_f .or. any(x /= x_f) .or. any(diffusion /= diffusion_f) .or. lambda /= lambda_f) stop 2 + if (any(h /= h_rt*gas_constant*t/molecular_weights)) stop 3 + write(*,'(*(ES25.16E3,1X))') cp_f, cv_f, e_f, lambda_f, diffusion_f, h + end do +end program +""" + + +@pytest.mark.parametrize("mechanism", MECHANISMS[:2]) +def test_fused_routines(tmp_path, mechanism): + """Caller-shaped routines equal the separate calls bitwise and agree with Cantera.""" + gas = ct.Solution(mechanism) + executable = compile_kernel(tmp_path, gas, driver_source=FUSED_DRIVER) + states = list(reference_states(gas)) + inputs = "\n".join(" ".join(map(str, [t, p, *y])) for t, p, y in states) + "\n" + result = subprocess.run([str(executable)], input=inputs, capture_output=True, text=True, check=True) + rows = np.array([np.fromstring(line, sep=" ") for line in result.stdout.splitlines()]) + assert len(rows) == len(states) + for actual, (t, p, y) in zip(rows, states): + gas.TPY = t, p, y + diffusion = gas.mix_diff_coeffs.copy() + for k in np.flatnonzero(y == 1): + diffusion[k] = gas.binary_diff_coeffs[k, k] + h = gas.standard_enthalpies_RT * ct.gas_constant * t / gas.molecular_weights + expected = np.concatenate(([gas.cp_mass, gas.cv_mass, gas.int_energy_mass, gas.thermal_conductivity], diffusion, h)) + np.testing.assert_allclose(actual, expected, rtol=2e-11, atol=1e-10) diff --git a/toolchain/mfc/thermochem/module.f90.mako b/toolchain/mfc/thermochem/module.f90.mako index 0a364e4f6..392366397 100644 --- a/toolchain/mfc/thermochem/module.f90.mako +++ b/toolchain/mfc/thermochem/module.f90.mako @@ -774,3 +774,131 @@ contains end subroutine get_species_mass_diffusivities_mixavg end module ${module_name} + +!> Routines shaped like MFC's call sites, kept out of the module above: +!> nearly every MFC file reads it, and nvfortran <= 24.7 crashes (fort1) +!> compiling unrelated files once it grows by these routines. +module ${module_name}_state + + use ${module_name} + + implicit none + + private + public :: get_mixture_caloric_state, get_species_enthalpies_mass, get_mixavg_transport_state + +contains + + !> Species cp/R and mixture cp, cv, e [J/kg] from one NASA7 pass. + !> A leaf routine: CCE faults on thermochem calls nested below a kernel. + subroutine get_mixture_caloric_state(temperature, mass_fractions, cp0_r, cp_mix, cv_mix, e_mix) + + GPU_ROUTINE(get_mixture_caloric_state) + + ${real_type}, intent(in) :: temperature + ${real_type}, intent(in), dimension(${sol.n_species}) :: mass_fractions + ${real_type}, intent(out), dimension(${sol.n_species}) :: cp0_r + ${real_type}, intent(out) :: cp_mix, cv_mix, e_mix + + ${real_type}, dimension(${sol.n_species}) :: shifted + + %for i, sp in enumerate(sol.species()): + cp0_r(${i+1}) = ${cgm(ce.poly_to_expr(sp.thermo, "temperature"))} + %endfor +${mass_average("cp_mix", "cp0_r")} + cp_mix = cp_mix * gas_constant + + shifted = cp0_r - 1.e0_${kind} +${mass_average("cv_mix", "shifted")} + cv_mix = cv_mix * gas_constant + + %for i, sp in enumerate(sol.species()): + shifted(${i+1}) = ${cgm(ce.poly_to_enthalpy_expr(sp.thermo, "temperature"))} + %endfor + shifted = shifted - 1.e0_${kind} +${mass_average("e_mix", "shifted")} + e_mix = e_mix * gas_constant * temperature + + end subroutine get_mixture_caloric_state + + !> Species enthalpies per unit mass [J/kg]. + subroutine get_species_enthalpies_mass(temperature, enthalpies) + + GPU_ROUTINE(get_species_enthalpies_mass) + + ${real_type}, intent(in) :: temperature + ${real_type}, intent(out), dimension(${sol.n_species}) :: enthalpies + + call get_species_enthalpies_rt(temperature, enthalpies) + %for i in range(sol.n_species): + enthalpies(${i+1}) = enthalpies(${i+1})*gas_constant*temperature/molecular_weights(${i+1}) + %endfor + + end subroutine get_species_enthalpies_mass + + !> Mixture-averaged transport of one state; shares the composition work. + !> The shared bodies are expanded, not called, to keep the call depth. + subroutine get_mixavg_transport_state(pressure, temperature, mass_fractions, & + mix_mol_weight, mole_fractions, mass_diffusivities_mixavg, & + mixture_thermal_conductivity_mixavg) + + GPU_ROUTINE(get_mixavg_transport_state) + + ${real_type}, intent(in) :: pressure, temperature + ${real_type}, intent(in), dimension(${sol.n_species}) :: mass_fractions + ${real_type}, intent(out) :: mix_mol_weight + ${real_type}, intent(out), dimension(${sol.n_species}) :: mole_fractions + ${real_type}, intent(out), dimension(${sol.n_species}) :: & + mass_diffusivities_mixavg + ${real_type}, intent(out) :: mixture_thermal_conductivity_mixavg + + ${real_type}, dimension(${sol.n_species}) :: conductivities, x_sum, denom + ${real_type}, dimension(${sol.n_species}, ${sol.n_species}) :: bdiff_ij + + call get_mixture_molecular_weight(mass_fractions, mix_mol_weight) + call get_mole_fractions(mix_mol_weight, mass_fractions, mole_fractions) +${diffusivities_body("mass_diffusivities_mixavg")} +${conductivity_body("mixture_thermal_conductivity_mixavg")} + end subroutine get_mixavg_transport_state + +end module ${module_name}_state +<%def name="mass_average(result, prop)">\ + ${result} = ( & + %for i in range(sol.n_species): + + inv_molecular_weights(${i+1})*mass_fractions(${i+1}) & + *${prop}(${i+1}) & + %endfor + ) +\ +<%def name="conductivity_body(result)">\ + call get_species_thermal_conductivities(temperature, conductivities) + + ${result} = 0.5_${kind}*(& + sum(mole_fractions*conductivities) + & + 1/sum(mole_fractions/conductivities)) +\ +<%def name="diffusivities_body(result)">\ + call get_species_binary_mass_diffusivities(temperature, bdiff_ij) + + %for sp in range(sol.n_species): + x_sum(${sp + 1}) = ${cgm(ce.diffusivity_mixture_rule_denom_expr( + sol, sp, Variable("mole_fractions"), Variable("bdiff_ij")))} + %endfor + + %for sp in range(sol.n_species): + denom(${sp + 1}) = x_sum(${sp + 1}) - & + mole_fractions(${sp + 1})/bdiff_ij(${sp + 1}, ${sp + 1}) + %endfor + + %for sp in range(sol.n_species): + if (denom(${sp + 1}) .gt. 0e0_${kind}) then + ${result}(${sp + 1}) = & + (mix_mol_weight - & + mole_fractions(${sp + 1})*molecular_weights(${sp + 1}))& + /(pressure * mix_mol_weight * denom(${sp + 1})) + else + ${result}(${sp + 1}) = & + bdiff_ij(${sp + 1}, ${sp + 1}) / pressure + end if + %endfor +\