From 28adc0d03e81e77afacc911c22a6a6df7178733d Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Tue, 23 Jan 2024 15:50:21 +0000 Subject: [PATCH] Add stretching for LMN --- src/navier.f90 | 44 +++++++++++++++++++++++++++++++++++++++----- src/transeq.f90 | 15 ++++++++++++++- 2 files changed, 53 insertions(+), 6 deletions(-) diff --git a/src/navier.f90 b/src/navier.f90 index 60d8ac8da..2f0b9b537 100644 --- a/src/navier.f90 +++ b/src/navier.f90 @@ -843,11 +843,11 @@ SUBROUTINE calc_divu_constraint(divu3, rho1, phi1) USE param, ONLY : ibirman_eos USE param, ONLY : xnu, prandtl USE param, ONLY : one - USE param, ONLY : iimplicit + USE param, ONLY : iimplicit, istret USE variables USE var, ONLY : ta1, tb1, tc1, td1, di1 - USE var, ONLY : phi2, ta2, tb2, tc2, td2, te2, di2 + USE var, ONLY : phi2, ta2, tb2, tc2, td2, te2, tf2, di2 USE var, ONLY : phi3, ta3, tb3, tc3, td3, rho3, di3 USE param, only : zero IMPLICIT NONE @@ -858,6 +858,8 @@ SUBROUTINE calc_divu_constraint(divu3, rho1, phi1) REAL(mytype), INTENT(IN), DIMENSION(xsize(1), xsize(2), xsize(3), numscalar) :: phi1 REAL(mytype), INTENT(OUT), DIMENSION(zsize(1), zsize(2), zsize(3)) :: divu3 + integer :: i, j, k + IF (ilmn.and.(.not.ibirman_eos)) THEN !!------------------------------------------------------------------------------ !! X-pencil @@ -901,6 +903,16 @@ SUBROUTINE calc_divu_constraint(divu3, rho1, phi1) tmp = iimplicit iimplicit = 0 CALL deryy (tc2, ta2, di2, sy, sfyp, ssyp, swyp, ysize(1), ysize(2), ysize(3), 1, zero) + if (istret /= 0) then + call dery (td2, ta2, di2, sy, ffyp, fsyp, fwyp, ppy, ysize(1), ysize(2), ysize(3), 1, zero) + do k = 1,ysize(3) + do j = 1,ysize(2) + do i = 1,ysize(1) + tc2(i,j,k) = tc2(i,j,k)*pp2y(j)-pp4y(j)*td2(i,j,k) + enddo + enddo + enddo + end if iimplicit = tmp IF (imultispecies) THEN tc2(:,:,:) = (xnu / prandtl) * tc2(:,:,:) / ta2(:,:,:) @@ -919,6 +931,16 @@ SUBROUTINE calc_divu_constraint(divu3, rho1, phi1) tmp = iimplicit iimplicit = 0 CALL deryy (td2, phi2(:,:,:,is), di2, sy, sfyp, ssyp, swyp, ysize(1), ysize(2), ysize(3), 1, zero) + if (istret /= 0) then + call dery (tf2, phi2(:,:,:,is), di2, sy, ffyp, fsyp, fwyp, ppy, ysize(1), ysize(2), ysize(3), 1, zero) + do k = 1,ysize(3) + do j = 1,ysize(2) + do i = 1,ysize(1) + td2(i,j,k) = td2(i,j,k)*pp2y(j)-pp4y(j)*tf2(i,j,k) + enddo + enddo + enddo + end if iimplicit = tmp tc2(:,:,:) = tc2(:,:,:) + (xnu / sc(is)) * (te2(:,:,:) / mol_weight(is)) * td2(:,:,:) ENDIF @@ -1041,13 +1063,13 @@ SUBROUTINE birman_drhodt_corr(drhodt1_next, rho1) USE decomp_2d, ONLY : mytype, xsize, ysize, zsize USE decomp_2d, ONLY : transpose_x_to_y, transpose_y_to_z, transpose_z_to_y, transpose_y_to_x - USE variables, ONLY : derxx, deryy, derzz + USE variables, ONLY : derxx, dery, deryy, derzz USE param, ONLY : nrhotime USE param, ONLY : xnu, prandtl - USE param, ONLY : iimplicit + USE param, ONLY : iimplicit, istret USE var, ONLY : td1, te1, di1, sx, sfxp, ssxp, swxp - USE var, ONLY : rho2, ta2, tb2, di2, sy, sfyp, ssyp, swyp + USE var, ONLY : rho2, ta2, tb2, tc2, di2, sy, sfyp, ssyp, swyp, ffyp, fsyp, fwyp, ppy, pp4y, pp2y USE var, ONLY : rho3, ta3, di3, sz, sfzp, sszp, swzp USE param, only : zero IMPLICIT NONE @@ -1057,6 +1079,8 @@ SUBROUTINE birman_drhodt_corr(drhodt1_next, rho1) REAL(mytype) :: invpe + integer :: i, j, k + invpe = xnu / prandtl CALL transpose_x_to_y(rho1(:,:,:,1), rho2) @@ -1068,6 +1092,16 @@ SUBROUTINE birman_drhodt_corr(drhodt1_next, rho1) iimplicit = -iimplicit CALL deryy (ta2,rho2,di2,sy,sfyp,ssyp,swyp,ysize(1),ysize(2),ysize(3),1, zero) + if (istret /= 0) then + call dery (tc2, rho2, di2, sy, ffyp, fsyp, fwyp, ppy, ysize(1), ysize(2), ysize(3), 1, zero) + do k = 1,ysize(3) + do j = 1,ysize(2) + do i = 1,ysize(1) + ta2(i,j,k) = ta2(i,j,k)*pp2y(j)-pp4y(j)*tc2(i,j,k) + enddo + enddo + enddo + end if iimplicit = -iimplicit ta2(:,:,:) = ta2(:,:,:) + tb2(:,:,:) CALL transpose_y_to_x(ta2, te1) diff --git a/src/transeq.f90 b/src/transeq.f90 index c392679aa..ee9f3ab98 100644 --- a/src/transeq.f90 +++ b/src/transeq.f90 @@ -1181,10 +1181,11 @@ subroutine continuity_rhs_eq(drho1, rho1, ux1, divu3) use param, only : ntime, nrhotime, ibirman_eos, zero use param, only : xnu, prandtl use param, only : iimplicit + use param, only : istret use variables use var, only : ta1, tb1, di1 - use var, only : rho2, uy2, ta2, tb2, di2 + use var, only : rho2, uy2, ta2, tb2, tc2, di2 use var, only : rho3, uz3, ta3, di3 implicit none @@ -1197,6 +1198,8 @@ subroutine continuity_rhs_eq(drho1, rho1, ux1, divu3) real(mytype) :: invpe + integer :: i, j, k + invpe = xnu / prandtl !! XXX All variables up to date - no need to transpose @@ -1219,6 +1222,16 @@ subroutine continuity_rhs_eq(drho1, rho1, ux1, divu3) iimplicit = -iimplicit call deryy (ta2,rho2,di2,sy,sfyp,ssyp,swyp,ysize(1),ysize(2),ysize(3),1, zero) + if (istret /= 0) then + call dery (tc2, rho2, di2, sy, ffyp, fsyp, fwyp, ppy, ysize(1), ysize(2), ysize(3), 1, zero) + do k = 1,ysize(3) + do j = 1,ysize(2) + do i = 1,ysize(1) + ta2(i,j,k) = ta2(i,j,k)*pp2y(j)-pp4y(j)*tc2(i,j,k) + enddo + enddo + enddo + end if iimplicit = -iimplicit ta2(:,:,:) = ta2(:,:,:) + tb2(:,:,:) call transpose_y_to_x(ta2, ta1)