Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add stretching for LMN #249

Draft
wants to merge 2 commits into
base: master
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
46 changes: 41 additions & 5 deletions src/navier.f90
Original file line number Diff line number Diff line change
Expand Up @@ -869,11 +869,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
Expand All @@ -884,6 +884,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
Expand Down Expand Up @@ -927,6 +929,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(:,:,:)
Expand All @@ -945,6 +957,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
Expand Down Expand Up @@ -1064,13 +1086,15 @@ SUBROUTINE extrapol_drhodt(drhodt1_next, rho1, drho1)

SUBROUTINE birman_drhodt_corr(drhodt1_next, rho1)

USE variables, ONLY : derxx, deryy, derzz
USE decomp_2d, ONLY : 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, 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
Expand All @@ -1080,6 +1104,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)
Expand All @@ -1091,6 +1117,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)
Expand Down
15 changes: 14 additions & 1 deletion src/transeq.f90
Original file line number Diff line number Diff line change
Expand Up @@ -1189,10 +1189,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
Expand All @@ -1205,6 +1206,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
Expand All @@ -1227,6 +1230,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)
Expand Down