Skip to content

forces: compute the unsteady contribution using the momentum-rhs #330

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.

Sign up for GitHub

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

Draft
mathrack wants to merge 2 commits into
base: master
Choose a base branch
from
Draft

forces: compute the unsteady contribution using the momentum-rhs #330

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
e26a31d
forces: compute the unsteady contribution using the momentum-rhs inst…
mathrack Dec 9, 2024
1192570
Fix force estimation
mathrack Dec 9, 2024
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
5 changes: 1 addition & 4 deletions src/case.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -388,10 +388,7 @@ subroutine postprocess_case(rho,ux,uy,uz,pp,phi,ep)

endif

if (iforces.eq.1) then
call force(ux,uy,uz,ep)
call restart_forces(1)
endif
if (iforces.eq.1) call force(ux,uy,uz,ep)

end subroutine postprocess_case
!##################################################################
Expand Down
125 changes: 29 additions & 96 deletions src/forces.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,37 +21,33 @@ module forces
implicit none

integer,save :: nvol,iforces,i2dsim
real(mytype),save,allocatable,dimension(:,:,:) :: ux01, uy01, ux11, uy11, ppi1
real(mytype),save,allocatable,dimension(:,:,:) :: dux, duy, ppi1
real(mytype),save,allocatable,dimension(:) :: xld,xrd,yld,yud,zfr,zbk
integer,save,allocatable,dimension(:) :: icvlf,icvrt,jcvlw,jcvup,kcvfr,kcvbk
integer,save,allocatable,dimension(:) :: icvlf_lx,icvrt_lx,icvlf_ly,icvrt_ly
integer,save,allocatable,dimension(:) :: jcvlw_lx,jcvup_lx,jcvlw_ly,jcvup_ly
integer,save,allocatable,dimension(:) :: kcvfr_lx,kcvbk_lx,kcvfr_ly,kcvbk_ly

character(len=*), parameter :: io_restart_forces = "restart-forces-io", &
resfile = "restart-forces"
private
public :: iforces, nvol, ppi1, init_forces, setup_forces, forces_unst, force

contains

subroutine init_forces

USE decomp_2d_io, only : decomp_2d_register_variable, decomp_2d_init_io
USE param
USE variables
implicit none

integer :: iv,stp1,stp2,h

call alloc_x(ux01)
call alloc_x(uy01)
call alloc_x(ux11)
call alloc_x(uy11)
call alloc_x(dux)
call alloc_x(duy)
call alloc_x(ppi1)

ux01 = zero
uy01 = zero
ux11 = zero
uy11 = zero
dux = zero
duy = zero
ppi1 = zero

allocate(icvlf(nvol),icvrt(nvol),jcvlw(nvol),jcvup(nvol),kcvfr(nvol),kcvbk(nvol))
allocate(icvlf_lx(nvol),icvrt_lx(nvol),icvlf_ly(nvol),icvrt_ly(nvol))
Expand Down Expand Up @@ -141,12 +137,6 @@ subroutine init_forces
endif
endif

call decomp_2d_init_io(io_restart_forces)
call decomp_2d_register_variable(io_restart_forces, "ux01", 1, 0, 0, mytype)
call decomp_2d_register_variable(io_restart_forces, "uy01", 1, 0, 0, mytype)
call decomp_2d_register_variable(io_restart_forces, "ux11", 1, 0, 0, mytype)
call decomp_2d_register_variable(io_restart_forces, "uy11", 1, 0, 0, mytype)

end subroutine init_forces

!
Expand Down Expand Up @@ -192,48 +182,26 @@ subroutine setup_forces(iounit)

end subroutine setup_forces

subroutine restart_forces(itest1)

USE decomp_2d_io
USE variables
USE param
USE MPI
! Store du/dt to compute the unsteady contribution later
subroutine forces_unst(rhs_dux, rhs_duy, px, py, dt)

implicit none
! Arguments
real(mytype), dimension(xsize(1),xsize(2),xsize(3)), intent(in) :: rhs_dux, rhs_duy, px, py
real(mytype), intent(in) :: dt

! Argument
integer, intent(in) :: itest1
! Local variables
integer :: i, j, k

! Exit if writing and invalid time step
if (itest1 == 1 .and. mod(itime, icheckpoint).ne.0) then
return
end if

if (itest1 == 1) then
call decomp_2d_open_io(io_restart_forces, resfile, decomp_2d_write_mode)
else
call decomp_2d_open_io(io_restart_forces, resfile, decomp_2d_read_mode)
endif
call decomp_2d_start_io(io_restart_forces, resfile)

if (itest1==1) then
!write
call decomp_2d_write_one(1,ux01,resfile,"ux01",0,io_restart_forces)
call decomp_2d_write_one(1,uy01,resfile,"uy01",0,io_restart_forces)
call decomp_2d_write_one(1,ux11,resfile,"ux11",0,io_restart_forces)
call decomp_2d_write_one(1,uy11,resfile,"uy11",0,io_restart_forces)
else
!read
call decomp_2d_read_one(1,ux01,resfile,"ux01",io_restart_forces)
call decomp_2d_read_one(1,uy01,resfile,"uy01",io_restart_forces)
call decomp_2d_read_one(1,ux11,resfile,"ux11",io_restart_forces)
call decomp_2d_read_one(1,uy11,resfile,"uy11",io_restart_forces)
endif
do k = 1, xsize(3)
do j = 1, xsize(2)
do i = 1, xsize(1)
dux(i,j,k) = dt*rhs_dux(i,j,k) - px(i,j,k)
duy(i,j,k) = dt*rhs_duy(i,j,k) - py(i,j,k)
end do
end do
end do

call decomp_2d_end_io(io_restart_forces, resfile)
call decomp_2d_close_io(io_restart_forces, resfile)

end subroutine restart_forces
end subroutine forces_unst

subroutine force(ux1,uy1,uz1,ep1)

Expand Down Expand Up @@ -287,28 +255,6 @@ subroutine force(ux1,uy1,uz1,ep1)
endif
enddo

if (itime.eq.1) then
do k = 1, xsize(3)
do j = 1, xsize(2)
do i = 1, xsize(1)
ux11(i,j,k)=ux1(i,j,k)
uy11(i,j,k)=uy1(i,j,k)
enddo
enddo
enddo
return
elseif (itime.eq.2) then
do k = 1, xsize(3)
do j = 1, xsize(2)
do i = 1, xsize(1)
ux01(i,j,k)=ux1(i,j,k)
uy01(i,j,k)=uy1(i,j,k)
enddo
enddo
enddo
return
endif

call derx (ta1,ux1,di1,sx,ffx,fsx,fwx,xsize(1),xsize(2),xsize(3),0,ubcx) ! dudx
call derx (tb1,uy1,di1,sx,ffxp,fsxp,fwxp,xsize(1),xsize(2),xsize(3),1,ubcy) ! dvdx
call transpose_x_to_y(ta1,ta2) ! dudx
Expand Down Expand Up @@ -370,12 +316,12 @@ subroutine force(ux1,uy1,uz1,ep1)
! relative to the volume, and, therefore, this has a sense
! of a "source".
! fac = (1.5*ux1(i,j,k)-2.0*ux01(i,j,k)+0.5*ux11(i,j,k))*epcv1(i,j,k)
fac = (onepfive*ux1(i,j,k)-two*ux01(i,j,k)+half*ux11(i,j,k))*(one-ep1(i,j,k))
fac = dux(i,j,k)*(one-ep1(i,j,k))
tsumx = tsumx+fac*dx*del_y(j+(xstart(2)-1))/dt !tsumx+fac*dx*dy/dt
!sumx(k) = sumx(k)+dudt1*dx*dy

! fac = (1.5*uy1(i,j,k)-2.0*uy01(i,j,k)+0.5*uy11(i,j,k))*epcv1(i,j,k)
fac = (onepfive*uy1(i,j,k)-two*uy01(i,j,k)+half*uy11(i,j,k))*(one-ep1(i,j,k))
fac = duy(i,j,k)*(one-ep1(i,j,k))
tsumy = tsumy+fac*dx*del_y(j+(xstart(2)-1))/dt !tsumy+fac*dx*dy/dt
!sumy(k) = sumy(k)+dudt1*dx*dy
enddo
Expand Down Expand Up @@ -634,12 +580,12 @@ subroutine force(ux1,uy1,uz1,ep1)
! relative to the volume, and, therefore, this has a sense
! of a "source".
! fac = (1.5*ux1(i,j,k)-2.0*ux01(i,j,k)+0.5*ux11(i,j,k))*epcv1(i,j,k)
fac = (onepfive*ux1(i,j,k)-two*ux01(i,j,k)+half*ux11(i,j,k))*(one-ep1(i,j,k))
fac = dux(i,j,k)*(one-ep1(i,j,k))
tsumx = tsumx+fac*dx*del_y(j+(xstart(2)-1))*dz/dt !tsumx+fac*dx*dy*dz/dt
!sumx(k) = sumx(k)+dudt1*dx*dy

! fac = (1.5*uy1(i,j,k)-2.0*uy01(i,j,k)+0.5*uy11(i,j,k))*epcv1(i,j,k)
fac = (onepfive*uy1(i,j,k)-two*uy01(i,j,k)+half*uy11(i,j,k))*(one-ep1(i,j,k))
fac = duy(i,j,k)*(one-ep1(i,j,k))
tsumy = tsumy+fac*dx*del_y(j+(xstart(2)-1))*dz/dt !tsumy+fac*dx*dy*dz/dt
!sumy(k) = sumy(k)+dudt1*dx*dy
enddo
Expand Down Expand Up @@ -920,19 +866,6 @@ subroutine force(ux1,uy1,uz1,ep1)
endif
endif

do k = 1, xsize(3)
do j = 1, xsize(2)
do i = 1, xsize(1)
ux11(i,j,k)=ux01(i,j,k)
uy11(i,j,k)=uy01(i,j,k)
ux01(i,j,k)=ux1(i,j,k)
uy01(i,j,k)=uy1(i,j,k)
enddo
enddo
enddo

return

end subroutine force

end module forces
9 changes: 3 additions & 6 deletions src/xcompact3d.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@ program xcompact3d
solve_poisson_mhd
use param, only : mhd_active
use particle, only : intt_particles
use forces, only : iforces, forces_unst

implicit none

Expand Down Expand Up @@ -56,6 +57,7 @@ program xcompact3d
endif
endif
call calculate_transeq_rhs(drho1,dux1,duy1,duz1,dphi1,rho1,ux1,uy1,uz1,ep1,phi1,divu3)
if (iforces.eq.1) call forces_unst(dux1,duy1,px1,py1,gdt(itr))

#ifdef DEBG
call check_transients()
Expand Down Expand Up @@ -216,12 +218,7 @@ subroutine init_xcompact3d()
call body(ux1,uy1,uz1,ep1)
endif

if (iforces.eq.1) then
call init_forces()
if (irestart==1) then
call restart_forces(0)
endif
endif
if (iforces.eq.1) call init_forces()

!####################################################################
! initialise mhd
Expand Down
Loading