module mod_ppm !! The **one-dimensional implementation** of the PPM !! (piecewise parabolic method) implicit none private #ifdef REAL_KIND integer, parameter :: fp = REAL_KIND #else integer, parameter :: fp = 8 #endif real(fp), parameter :: eps = epsilon(1.0_fp) !! 极小值 public :: fp, cal_cfl_time_step, adv1d_by_ppm contains subroutine cal_cfl_time_step(dx, u, dt, dt_, nt) !! 计算满足 CFL(Courant–Friedrichs–Lewy) 条件的时间积分步长 real(fp), intent(in) :: dx(:) !! 网格分辨率: m real(fp), intent(in) :: u(:) !! 风速: m/s real(fp), intent(in) :: dt !! 全局迭代步长: s real(fp), intent(out) :: dt_ !! 局部迭代步长: s integer, intent(out) :: nt !! 迭代多少次 ! local variable real(fp), parameter :: factor = 0.5_fp !! 满足CFL条件的程度 < 0.5 ! 计算 CFL 约束的最小时间步长 dt_ = factor * minval(dx/abs(u)) ! 确保 dt_ 能整除 dt，并计算迭代次数 if (dt_ < dt) then nt = int(dt/dt_) + 1 dt_ = dt/nt else dt_ = dt nt = 1 end if end subroutine cal_cfl_time_step subroutine adv1d_by_ppm(dt, u, c, advc, dx, ds, volume) !! 基于PPM方法实现的一维平流函数 ! input args real(fp), intent(in) :: dt !! 时间间隔: s real(fp), intent(in) :: u(:) !! 网格右边界风速: m/s ! output args real(fp), intent(inout) :: c(:) !! 网格浓度: ug/m3 real(fp), intent(out) :: advc(:) !! 浓度变化: ug/m3 ! optional args: dx 和 ds必须有一个 real(fp), optional, intent(in) :: dx !! 网格分辨率: m real(fp), optional, intent(in) :: ds(:) !! 网格分辨率: m real(fp), optional, intent(in) :: volume(:) !! 体积校正因子，dx*dy*dz; ! local variables ! mass grid real(fp) :: c6(size(u)) !! 决定了抛物线的凸凹特征 real(fp) :: deltaC(size(u)) !! rightC - leftC ! stag grid: 第i个mass网格的右边为stag网格的i real(fp) :: egdeC(0:size(u)) !! 格子的边界处的浓度: ug/m3 real(fp) :: leftC(size(u)) !! 网格的左边界处的浓度 leftC = (0:n-1): ug/m3 real(fp) :: rightC(size(u)) !! 网格的右边界处的浓度 rightC = (1:n): ug/m3 real(fp) :: fluidC(size(u)) !! 边界处的浓度随着风场输出(一小段距离)的平均浓度: ug/m3 real(fp) :: flux(0:size(u)) !! 边界处的通量: ug/m^2/s real(fp) :: CFL !! Courant-Friedrichs-Lewy条件（u*dt/dx） real(fp) :: fL !! 网格的左边界处通量: ug/m^2/s real(fp) :: fR !! 网格的右边界处通量: ug/m^2/s integer :: i, n n = size(u) if (present(dx)) then call cal_conc_at_uniform_mesh_edge(n, c, egdeC) else if (present(ds)) then call cal_ppm_non_uniform_mesh_edge_conc(n, ds, c, egdeC) else write(*,*) 'You must give a dx' stop 1 end if leftC = egdeC(0:n-1) rightC = egdeC(1:n) ! step 3: 计算抛物线参数，并对边界处浓度进行约束 do i=1, n ! 更新边界处浓度 if ((rightC(i) - c(i))*(c(i) - leftC(i)) > 0.) then deltaC(i) = rightC(i) - leftC(i) c6(i) = 6.*(c(i) - (rightC(i) + leftC(i))/2.) ! Overshoot cases if (deltaC(i)*c6(i) > deltaC(i)*deltaC(i)) then leftC(i) = 3.*c(i) - 2.*rightC(i) else if (-deltaC(i)*deltaC(i) > deltaC(i)*c6(i)) then rightC(i) = 3.*c(i) - 2.*leftC(i) endif else ! 局地极值 leftC(i) = c(i) rightC(i) = c(i) end if ! 确定每个格子抛物线参数：rightC, leftC, c6; deltaC(i) = rightC(i) - leftC(i) c6(i) = 6.*(c(i) - (rightC(i) + leftC(i))/2.) end do ! step 4: Compute fluxes ! fluxes from boundary cells flux(0) = 0 ! surface boundary: impermeable boundary condition if (u(1) >= 0.) flux(1) = u(1) * rightC(1) ! First order if (u(1) < 0.) flux(1) = u(1) * leftC(1+1) ! First order if (u(n-1) >= 0.) flux(n-1) = u(n-1) * rightC(n-1) ! First order if (u(n-1) < 0.) flux(n-1) = u(n-1) * leftC(n) ! First order if (u(n) > 0.) flux(n) = u(n) * rightC(n) ! 顶边界, 只出去不进 ! fluxes from the parabolic distribution CFL = 0.5 ! 初始化，没有用 do i = 2, n-2 if (present(dx)) then CFL = abs((u(i)*dt)/dx) ! Guarantee y = u(i)*dt > 0 else if (present(ds)) then CFL = abs((u(i)*dt)/ds(i)) ! Guarantee y = u(i)*dt > 0 ! 这里的ds, 是不是应该是 (ds(i) + ds(i))/2 end if if (u(i) > 0.) then ! 公式 2.12 fluidC(i) = rightC(i) - CFL/2 * (deltaC(i) - (1. - CFL*2./3.)*c6(i)) else fluidC(i) = leftC(i+1) + CFL/2 * (deltaC(i+1) + (1. - CFL*2./3.)*c6(i+1)) end if flux(i) = u(i) * fluidC(i) ! 公式2.14 end do ! step 5: 时间积分, 前向欧拉 do i = 2, n-1 fR = flux(i) fL = flux(i-1) if (present(volume)) then ! 进行校正体积校正，保证质量守恒 if (fL > 0 .and. i>1) fL = fL*volume(i-1)/volume(i) if (fR < 0 .and. i<n) fR = fR*volume(i+1)/volume(i) end if if (present(dx)) then advc(i) = dt/dx * (fL - fR) else if (present(ds)) then advc(i) = dt/ds(i) * (fL - fR) end if c(i) = c(i) + advc(i) end do ! 极小值，数值误差 where(c < eps) c = 0. end subroutine adv1d_by_ppm subroutine cal_conc_at_uniform_mesh_edge(n, c, egdeC) !! 计算均匀网格边界处的浓度(二次抛物线) ! input args integer, intent(in) :: n !! 网格数 real(fp), intent(in) :: c(n) !! 网格浓度: ug/m3 ! output args real(fp), intent(out) :: egdeC(0:n) !! conc at mesh egde: ug/m3 ! local variables integer :: I real(fp) :: slopC(n) !! average slop ! step 1: 计算 2~n-1 个网格的平均坡度 (slopC) do i = 2, n-1 ! Compute average slope in the i'th cell slopC(i) = 0.5 * (c(i+1) - c(i-1)) ! Constraint: 确保网格边界处浓度在左右格子平均浓度之间 if ((c(i+1) - c(i))*(c(i) - c(i-1)) > 0.) then slopC(i) = sign(1._fp, slopC(i))* & min(abs(slopC(i)), & 2.*abs(c(i+1) - c(i)), & 2.*abs(c(i) - c(i-1))) else ! 局地极值情况 slopC(i) = 0 end if end do ! step 2: Compute concentrations at cell boundaries ! Zero order polynomial egdeC(0) = c(1) egdeC(n) = c(n) ! First order polynomial egdeC(1) = (c(2) + c(1))/2. egdeC(n-1) = (c(n) + c(n-1))/2. ! 求边界处的值 do i = 2, n-2 ! 中间网格 egdeC(i) = c(i) + 0.5*(c(i+1) - c(i)) + (slopC(i) - slopC(i+1))/6. end do end subroutine cal_conc_at_uniform_mesh_edge subroutine cal_ppm_non_uniform_mesh_edge_conc(n, dx, c, egdeC) !! 计算变网格边界处浓度((二次抛物线)) ! input args integer, intent(in) :: n !! 网格数 real(fp), intent(in) :: dx(n) !! 分辨率: m real(fp), intent(in) :: c(n) !! 网格浓度: ug/m3 ! output args real(fp), intent(out) :: egdeC(0:n) !! conc at mesh egde: ug/m3 ! local variables integer :: I real(fp) :: x0, x1, x2, x3 ! temp lattice vars. real(fp) :: alpha (2:n-1) ! temp lattice var. real(fp) :: chi (2:n-1) ! lattice var. for dc real(fp) :: psi (2:n-1) ! lattice var. for dc real(fp) :: mu (2:n-1) ! lattice var. for cm real(fp) :: nu (2:n-1) ! lattice var. for cm real(fp) :: lambda(2:n-1) ! lattice var. for cm real(fp) :: slopC(n) ! step 1: args do i = 2, n-1 alpha(i) = dx(i) + dx(i+1) x0 = dx(i-1) + dx(i) x2 = dx(i) / (x0 + dx(i+1)) chi(i) = x2 * (dx(i-1) + x0) / alpha(i) psi(i) = x2 * (alpha(i) + dx(i+1)) / x0 end do do i = 2, n-2 x1 = dx(i)/alpha(i) x2 = 2.0_fp * dx(i+1) / alpha(i) x3 = 1.0_fp / (dx(i-1) + alpha(i) + dx(i+2)) mu(i) = x3 * dx(i) * (dx(i-1) + dx(i)) / (dx(i) + alpha(i)) nu(i) = x3 * dx(i+1) * (dx(i+1) + dx(i+2)) / (dx(i+1) + alpha(i)) lambda(i) = x1 + mu(i) * x2 - 2.0*nu(i) * x1 end do ! step 1: 平均坡度 (slopC) do i = 2, n-1 ! Compute average slope in the i'th cell slopC(i) = chi(i)*(c(i+1) - c(i)) + psi(i)*(c(i) - c(i-1)) ! equation (1.7) ! Constraint: 确保梯度不要太大 if ((c(i+1) - c(i))*(c(i) - c(i-1)) > 0.) then ! equation (1.8) slopC(i) = sign(1._fp, slopC(i))* & min( abs(slopC(i)), & 2.*abs(c(i+1) - c(i)), & 2.*abs(c(i) - c(i-1)) ) else ! 局地极值情况 slopC(i) = 0 end if end do ! step 2: Compute concentrations at cell boundaries egdeC(0) = c(1) ! Zero order polynomial egdeC(n) = c(n) egdeC(1) = (dx(1)*c(2) + dx(2)*c(1)) / (dx(1) + dx(2)) ! First order polynomial egdeC(n-1) = (dx(n-1)*c(n) + dx(n)*c(n-1)) / (dx(n-1) + dx(n)) ! second order polynomial inside the domain do i = 2, n-1 ! equation (1.6) egdeC(i) = c(i) + lambda(i)*(c(i+1) - c(i)) - mu(i)*slopC(i+1) + nu(i)*slopC(i) end do end subroutine cal_ppm_non_uniform_mesh_edge_conc end module mod_ppm