594 lines
19 KiB
Fortran
594 lines
19 KiB
Fortran
subroutine rhs_scalars
|
|
|
|
use m_openmpi
|
|
use m_io
|
|
use m_parameters
|
|
use m_fields
|
|
use m_work
|
|
use x_fftw
|
|
use m_timing
|
|
use m_les
|
|
|
|
implicit none
|
|
|
|
integer :: i, j, k, n, n1, n2, nv, ns_lo, ns_hi
|
|
real*8 :: rtmp1, rtmp2, wnum2, r11, r12, r21, r22, r31, r32
|
|
|
|
! calculate turbulent viscosity, if any
|
|
if (les) call les_get_turb_visc
|
|
|
|
! If we're not advancing scalars, put the real-space velocities
|
|
! in wrk1...3 and return
|
|
! same is done if dealias=0, that is, we use 2/3 rule.
|
|
! if dealias=1 (phase shifts) then this is done later in the subroutine
|
|
|
|
! also if doing LES and n_les (the # of les-related auxilary scalars) is
|
|
! greater than 0, then we need to transport these LES-related scalars,
|
|
! even if n_scalars=0.
|
|
|
|
! thus we do the obligatory part (tranfer velocities to X-space) and quit
|
|
! only when we are not transporting any scalars and if there are no
|
|
! LES-related scalars.
|
|
if (.not.int_scalars .and. n_les .eq. 0) then
|
|
! converting velocities to the real space and returning
|
|
wrk(:,:,:,1:3) = fields(:,:,:,1:3)
|
|
do n = 1,3
|
|
call xFFT3d(-1,n)
|
|
end do
|
|
return
|
|
end if
|
|
|
|
! making the RHS for all scalars zero
|
|
wrk(:,:,:,4:3+n_scalars+n_les) = zip
|
|
|
|
!--------------------------------------------------------------------------------
|
|
! If dealias=0, performing the 2/3 rule dealiasing on scalars
|
|
!--------------------------------------------------------------------------------
|
|
|
|
if (dealias.eq.0) then
|
|
|
|
! converting velocities to the real space
|
|
wrk(:,:,:,1:3) = fields(:,:,:,1:3)
|
|
do n = 1,3
|
|
call xFFT3d(-1,n)
|
|
end do
|
|
|
|
! Do each scalar one at a time. Keep the velocities in wrk1:3 intact
|
|
! because they are needed later.
|
|
|
|
! first we need to know which scalars do we want to transport.
|
|
! There are three cases:
|
|
! (0) No LES extra scalars, and passive scalars have not been initialized yet
|
|
! Thus this subroutine calculates the IFFT of velocities and exits.
|
|
! This is taken care of earlier.
|
|
! (1) Both passive scalars and LES extra scalars are transported
|
|
! This is possible when n_les > 0 and int_scalars=.true.
|
|
! (2) Only LES extra scalars are transported
|
|
! This is possible if and only if (.not.int_scalars .and. n_les>0)
|
|
!
|
|
! The last two cases are taken care of by prescribing ns_lo and ns_hi, the
|
|
! smallest and largest number of the scalar that needs to be transported.
|
|
|
|
ns_lo = 1;
|
|
ns_hi = n_scalars + n_les
|
|
if (.not.int_scalars) ns_lo = n_scalars + 1
|
|
|
|
do n = ns_lo, ns_hi
|
|
|
|
wrk(:,:,:,0) = fields(:,:,:,3+n)
|
|
call xFFT3d(-1,0)
|
|
|
|
! Products of the scalar and velocities
|
|
do i = 1,3
|
|
wrk(:,:,:,n+2+i) = wrk(:,:,:,0) * wrk(:,:,:,i)
|
|
call xFFT3d(1,n+2+i)
|
|
end do
|
|
|
|
! Assembling the RHS in wrk(:,:,:,3+n)
|
|
do k = 1,nz
|
|
do j = 1,ny
|
|
do i = 1,nx+1,2
|
|
|
|
! If the dealiasing option is 2/3-rule (dealias=0) then we retain the modes
|
|
! inside the cube described by $| k_i | \leq k_{max}$, $i=1,2,3$.
|
|
! The rest of the modes is purged
|
|
|
|
if (ialias(i,j,k) .gt. 0) then
|
|
! all the wavenumbers that are greater than kmax get zeroed out
|
|
wrk(i ,j,k,3+n) = zip
|
|
wrk(i+1,j,k,3+n) = zip
|
|
|
|
else
|
|
! taking the convective term, multiply it by "i"
|
|
! (see how it's done in x_fftw.f90)
|
|
! and adding the diffusion term
|
|
|
|
! also using the fact that the waveunmbers for (i,j,k) are the same
|
|
! as wavenumbers for (i+1,j,k)
|
|
|
|
! i * (a + ib) + d = -b + ia + d
|
|
rtmp1 = akx(i+1)*wrk(i+1,j,k,3+n) + aky(k)*wrk(i+1,j,k,4+n) + akz(j)*wrk(i+1,j,k,5+n)
|
|
rtmp2 = akx(i )*wrk(i ,j,k,3+n) + aky(k)*wrk(i ,j,k,4+n) + akz(j)*wrk(i ,j,k,5+n)
|
|
|
|
wnum2 = akx(i)**2 + aky(k)**2 + akz(j)**2
|
|
|
|
wrk(i ,j,k,3+n) = rtmp1 - pe(n) * wnum2*fields(i ,j,k,3+n)
|
|
wrk(i+1,j,k,3+n) = - rtmp2 - pe(n) * wnum2*fields(i+1,j,k,3+n)
|
|
|
|
end if
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
! Now adding the reaction part (for scalars only, not for LES-related quantities
|
|
! that are formally scalars with indicies n_scalars+1...n_scalars+n_les )
|
|
if (n .le. n_scalars) then
|
|
if (scalar_type(n).ge.100) then
|
|
call add_reaction(n)
|
|
call dealias_rhs(3+n)
|
|
end if
|
|
end if
|
|
|
|
end do
|
|
|
|
end if
|
|
|
|
!--------------------------------------------------------------------------------
|
|
! If dealias=1, performing the phase shift and truncation on scalars
|
|
! The main ideology is as follows. First we evaluate the phase-shifted
|
|
! quantities, then not phase shifted. This is done in order to have more
|
|
! working space: at the beginning we have all work arrays available, but at the
|
|
! end we must put sin/cos factors in wrk0 and u.v.w in real space in wrk1...3.
|
|
! They are going to be used again in rhs_velocity later.
|
|
!--------------------------------------------------------------------------------
|
|
phase_shifting_dealiasing: if (dealias.eq.1) then
|
|
|
|
! define the sin/cos factors that are used in phase shifting.
|
|
! computing sines and cosines for the phase shift of dx/2,dy/2,dz/2
|
|
! and putting them into wrk0
|
|
do k = 1,nz
|
|
do j = 1,ny
|
|
do i = 1,nx+1,2
|
|
wrk(i ,j,k,0) = cos(half*(akx(i )+aky(k)+akz(j))*dx)
|
|
wrk(i+1,j,k,0) = sin(half*(akx(i+1)+aky(k)+akz(j))*dx)
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
! Doing all scalars at the same time. We can do this because the
|
|
! number of work arrays that are available to us is 3+n+2. First
|
|
! three later will be taken by velocities, and the last two are
|
|
! primary work arrays here.
|
|
|
|
! First, get phase shifted velocities and scalars
|
|
do n = 1, 3 + n_scalars + n_les
|
|
! phase-shifting the quantity
|
|
do k = 1,nz
|
|
do j = 1,ny
|
|
do i = 1,nx+1,2
|
|
wrk(i ,j,k,n) = fields(i ,j,k,n) * wrk(i,j,k,0) - fields(i+1,j,k,n) * wrk(i+1,j,k,0)
|
|
wrk(i+1,j,k,n) = fields(i+1,j,k,n) * wrk(i,j,k,0) + fields(i ,j,k,n) * wrk(i+1,j,k,0)
|
|
end do
|
|
end do
|
|
end do
|
|
! transforming it to real space
|
|
call xFFT3d(-1,n)
|
|
end do
|
|
|
|
! now we have two vacant arrays: n_scalars+n_les+4 and n_scalars+n_les+5. Work in them
|
|
n1 = n_scalars + n_les + 4
|
|
n2 = n_scalars + n_les + 5
|
|
|
|
! do one scalar at a time
|
|
phase_shifted_rhs: do n = 4, 3 + n_scalars + n_les
|
|
|
|
! getting all three products of phase-shifted scalar and phase-shifted velocities
|
|
! using three work arrays: n, n1 and n2
|
|
wrk(:,:,:,n1) = wrk(:,:,:,n) * wrk(:,:,:,1)
|
|
wrk(:,:,:,n2) = wrk(:,:,:,n) * wrk(:,:,:,2)
|
|
wrk(:,:,:,n ) = wrk(:,:,:,n) * wrk(:,:,:,3)
|
|
! transforming them to Fourier space
|
|
call xFFT3d(1,n1)
|
|
call xFFT3d(1,n2)
|
|
call xFFT3d(1,n )
|
|
|
|
! phase shifting them back using wrk0 and adding -0.5*(ik) to the RHS for the scalar
|
|
do k = 1,nz
|
|
do j = 1,ny
|
|
do i = 1,nx+1,2
|
|
|
|
if (ialias(i,j,k) .gt. 1) then
|
|
wrk(i:i+1,j,k,n) = zip
|
|
else
|
|
! (u+)*(phi+) phase shifted back
|
|
r11 = wrk(i ,j,k,n1) * wrk(i,j,k,0) + wrk(i+1,j,k,n1) * wrk(i+1,j,k,0)
|
|
r12 = wrk(i+1,j,k,n1) * wrk(i,j,k,0) - wrk(i ,j,k,n1) * wrk(i+1,j,k,0)
|
|
! (v+)*(phi+) phase shifted back
|
|
r21 = wrk(i ,j,k,n2) * wrk(i,j,k,0) + wrk(i+1,j,k,n2) * wrk(i+1,j,k,0)
|
|
r22 = wrk(i+1,j,k,n2) * wrk(i,j,k,0) - wrk(i ,j,k,n2) * wrk(i+1,j,k,0)
|
|
! (w+)*(phi+) phase shifted back
|
|
r31 = wrk(i ,j,k,n ) * wrk(i,j,k,0) + wrk(i+1,j,k,n ) * wrk(i+1,j,k,0)
|
|
r32 = wrk(i+1,j,k,n ) * wrk(i,j,k,0) - wrk(i ,j,k,n ) * wrk(i+1,j,k,0)
|
|
! adding -0.5*(ik)*(the result) to the RHSs for the scalar
|
|
wrk(i ,j,k,n) = + 0.5d0 * ( akx(i+1)*r12 + aky(k)*r22 + akz(j)*r32 )
|
|
wrk(i+1,j,k,n) = - 0.5d0 * ( akx(i )*r11 + aky(k)*r21 + akz(j)*r31 )
|
|
end if
|
|
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
end do phase_shifted_rhs
|
|
|
|
! at this moment wrk4...3+n_scalars+n_les contain the half of the convective term, which was
|
|
! obtained from the phase shifted quantities. Now we need to add the other half of the
|
|
! convective term (the one that is obtained by multiplication of no-phase-shifted stuff)
|
|
|
|
! first get the velocities into the real space and put them in wrk1...3
|
|
! this should remain in there untouched, to be used in rhs_velocity later
|
|
do n = 1,3
|
|
wrk(:,:,:,n) = fields(:,:,:,n)
|
|
call xFFT3d(-1,n)
|
|
end do
|
|
|
|
! now do scalars one at a time since we don't have enough storage to do them
|
|
! all at once
|
|
|
|
not_phase_shifted_rhs: do n = 4, 3 + n_scalars + n_les
|
|
|
|
! get the scalar into the real space and put it in the wrk(0)
|
|
wrk(:,:,:,0) = fields(:,:,:,n)
|
|
call xFFT3d(-1,0)
|
|
|
|
! calculate the product of the scalar with velocitiy components
|
|
wrk(:,:,:,n1) = wrk(:,:,:,0) * wrk(:,:,:,1)
|
|
wrk(:,:,:,n2) = wrk(:,:,:,0) * wrk(:,:,:,2)
|
|
wrk(:,:,:, 0) = wrk(:,:,:,0) * wrk(:,:,:,3)
|
|
! transform it to the Fourier space
|
|
call xFFT3d(1,n1)
|
|
call xFFT3d(1,n2)
|
|
call xFFT3d(1, 0)
|
|
|
|
! add the -0.5*(ik)*results to the RHS, along with the diffusion term
|
|
do k = 1,nz
|
|
do j = 1,ny
|
|
do i = 1,nx+1,2
|
|
|
|
if (ialias(i,j,k) .lt. 2) then
|
|
|
|
rtmp1 = akx(i+1)*wrk(i+1,j,k,n1) + aky(k)*wrk(i+1,j,k,n2) + akz(j)*wrk(i+1,j,k,0)
|
|
rtmp2 = akx(i )*wrk(i ,j,k,n1) + aky(k)*wrk(i ,j,k,n2) + akz(j)*wrk(i ,j,k,0)
|
|
|
|
wnum2 = akx(i)**2 + aky(k)**2 + akz(j)**2
|
|
wrk(i ,j,k,n) = wrk(i ,j,k,n) + 0.5d0 * rtmp1 - pe(n-3) * wnum2*fields(i ,j,k,n)
|
|
wrk(i+1,j,k,n) = wrk(i+1,j,k,n) - 0.5d0 * rtmp2 - pe(n-3) * wnum2*fields(i+1,j,k,n)
|
|
|
|
end if
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
end do not_phase_shifted_rhs
|
|
|
|
! ------------------------------------------------------------
|
|
! Add the reaction rates to the RHS
|
|
!
|
|
! The LES-related scalars are not affected because the
|
|
! upper bound for n is 3+n_scalars, not 3+n_scalars+n_les
|
|
! -----------------------------------------------------------
|
|
|
|
reaction_rates: do n = 4, 3+n_scalars
|
|
|
|
if (scalar_type(n-3) .gt. 300) then
|
|
|
|
!!$ ! putting phase shifted scalar in wrk(n1)
|
|
!!$ do k = 1,nz
|
|
!!$ do j = 1,ny
|
|
!!$ do i = 1,nx+1,2
|
|
!!$ wrk(i ,j,k,n1) = fields(i ,j,k,n) * wrk(i,j,k,0) - fields(i+1,j,k,n) * wrk(i+1,j,k,0)
|
|
!!$ wrk(i+1,j,k,n1) = fields(i+1,j,k,n) * wrk(i,j,k,0) + fields(i ,j,k,n) * wrk(i+1,j,k,0)
|
|
!!$ end do
|
|
!!$ end do
|
|
!!$ end do
|
|
!!$ ! transforming it to real space
|
|
!!$ call xFFT3d(-1,n1)
|
|
!!$ ! putting the phase-shifted reaction rate in wrk(n2)
|
|
!!$ call scalar_reaction_rate(n1,n2)
|
|
!!$ ! transforming to Fourier space
|
|
!!$ call xFFT3d(1,n2)
|
|
!!$ ! phase shifting back and adding a half of it to the RHS
|
|
!!$ do k = 1,nz
|
|
!!$ do j = 1,ny
|
|
!!$ do i = 1,nx+1,2
|
|
!!$ if (ialias(i,j,k) .le. 1) then
|
|
!!$ ! phase shifting back
|
|
!!$ r11 = wrk(i ,j,k,n2) * wrk(i,j,k,0) + wrk(i+1,j,k,n2) * wrk(i+1,j,k,0)
|
|
!!$ r12 = wrk(i+1,j,k,n2) * wrk(i,j,k,0) - wrk(i ,j,k,n2) * wrk(i+1,j,k,0)
|
|
!!$ ! adding 0.5*(the result) to the RHSs for the scalar
|
|
!!$ wrk(i ,j,k,n) = wrk(i ,j,k,n) + 0.5d0 * r11
|
|
!!$ wrk(i+1,j,k,n) = wrk(i+1,j,k,n) + 0.5d0 * r12
|
|
!!$ end if
|
|
!!$ end do
|
|
!!$ end do
|
|
!!$ end do
|
|
!!$
|
|
!!$ ! scond part: doing the same thing with not-phase-shifted scalar
|
|
!!$ ! putting it in wrk(n1)
|
|
!!$ wrk(:,:,:,n1) = fields(:,:,:,n)
|
|
!!$ ! transforming it to real space
|
|
!!$ call xFFT3d(-1,n1)
|
|
!!$ ! putting the phase-shifted reaction rate in wrk(n2)
|
|
!!$ call scalar_reaction_rate(n1,n2)
|
|
!!$ ! transforming to Fourier space
|
|
!!$ call xFFT3d(1,n2)
|
|
!!$ ! phase shifting back and adding a half of it to the RHS
|
|
!!$ do k = 1,nz
|
|
!!$ do j = 1,ny
|
|
!!$ do i = 1,nx+1,2
|
|
!!$ if (ialias(i,j,k) .le. 1) then
|
|
!!$ ! phase shifting back
|
|
!!$ r11 = wrk(i ,j,k,n2) * wrk(i,j,k,0) + wrk(i+1,j,k,n2) * wrk(i+1,j,k,0)
|
|
!!$ r12 = wrk(i+1,j,k,n2) * wrk(i,j,k,0) - wrk(i ,j,k,n2) * wrk(i+1,j,k,0)
|
|
!!$ ! adding 0.5*(the result) to the RHSs for the scalar
|
|
!!$ wrk(i ,j,k,n) = wrk(i ,j,k,n) + 0.5d0 * r11
|
|
!!$ wrk(i+1,j,k,n) = wrk(i+1,j,k,n) + 0.5d0 * r12
|
|
!!$ end if
|
|
!!$ end do
|
|
!!$ end do
|
|
!!$ end do
|
|
|
|
if (scalar_type(n-3) .eq. 311) then
|
|
|
|
! since the reaction rate is cubic, we need to apply some severe truncation.
|
|
|
|
! putting the scalar in wrk(n1)
|
|
wrk(:,:,:,n1) = fields(:,:,:,n)
|
|
! truncating it so only the modes with |k_i| < nx/4 remain
|
|
do k = 1,nz
|
|
do j = 1,ny
|
|
do i = 1,nx+1,2
|
|
if (abs(akx(i)).gt.nx/4 .or. abs(aky(k)).gt.nx/4 .or. abs(akz(j)).gt.nx/4) then
|
|
wrk(i ,j,k,n1) = zip
|
|
wrk(i+1,j,k,n1) = zip
|
|
end if
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
! transforming it to real space
|
|
call xFFT3d(-1,n1)
|
|
|
|
! self-adjusting bistable reaction needs the mean value of the scalar
|
|
rtmp1 = fields(1,1,1,n)/nxyz_all
|
|
call MPI_BCAST(rtmp1, 1, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)
|
|
|
|
! getting the reaction rate
|
|
wrk(:,:,:,n2) = reac_sc(n-3) * (one - wrk(:,:,:,n1)**2) * (wrk(:,:,:,n1) - rtmp1)
|
|
|
|
! transforming it to Fourier space
|
|
call xFFT3d(1,n2)
|
|
|
|
! adding to the RHS
|
|
wrk(:,:,:,n) = wrk(:,:,:,n) + wrk(:,:,:,n2)
|
|
|
|
else
|
|
write(out,*) "The scalar type has a reaction rate that is not supported yet:", scalar_type(n-3)
|
|
call flush(out)
|
|
call my_exit(-1)
|
|
end if
|
|
end if
|
|
end do reaction_rates
|
|
|
|
end if phase_shifting_dealiasing
|
|
|
|
! special case - passive scalar with the uniform gradient as a source
|
|
! adding the source term - the first component of velocity, because we assume
|
|
! that the uniform gradient has slope 1 and direction in the x-direction
|
|
gradient_source: do n = 1, n_scalars
|
|
if (scalar_type(n) .eq. 0) wrk(:,:,:,n+3) = wrk(:,:,:,n+3) - fields(:,:,:,1)
|
|
end do gradient_source
|
|
|
|
|
|
|
|
!--------------------------------------------------------------------------------
|
|
! Add LES to the RHS of all the scalars
|
|
!--------------------------------------------------------------------------------
|
|
les_active: if (les) then
|
|
call les_rhs_scalars
|
|
end if les_active
|
|
|
|
|
|
return
|
|
end subroutine rhs_scalars
|
|
|
|
!================================================================================
|
|
!================================================================================
|
|
subroutine add_reaction(n)
|
|
|
|
use m_openmpi
|
|
use m_io
|
|
use m_parameters
|
|
use m_fields
|
|
use m_work
|
|
use x_fftw
|
|
|
|
implicit none
|
|
|
|
integer :: n, rtype
|
|
real*8 :: scmean, rrate
|
|
|
|
! reaction type
|
|
rtype = scalar_type(n)/100
|
|
|
|
! raction rate
|
|
rrate = reac_sc(n)
|
|
|
|
select case (rtype)
|
|
case (1)
|
|
|
|
! KPP reaction rate
|
|
wrk(:,:,:,0) = rrate * (1.d0 - wrk(:,:,:,0)**2)
|
|
|
|
case (2)
|
|
|
|
! symmetric bistable
|
|
wrk(:,:,:,0) = rrate * (1.d0 - wrk(:,:,:,0)**2) * wrk(:,:,:,0)
|
|
|
|
case (3)
|
|
|
|
! self-adjusting bistable
|
|
scmean = fields(1,1,1,3+n)/nxyz_all
|
|
call MPI_BCAST(scmean, 1, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)
|
|
|
|
wrk(:,:,:,0) = rrate * (1.d0 - wrk(:,:,:,0)**2) * &
|
|
(wrk(:,:,:,0) - scmean)
|
|
case default
|
|
|
|
write(out,*) "Unknown reaction rate"
|
|
call flush(out)
|
|
|
|
stop
|
|
|
|
end select
|
|
|
|
|
|
! FFT the reaction into the Fourier space
|
|
call xFFT3d(1,0)
|
|
|
|
! Adding reaction to the RHS in wrk(:,:,:,3+n)
|
|
|
|
wrk(:,:,:,3+n) = wrk(:,:,:,3+n) + wrk(:,:,:,0)
|
|
|
|
end subroutine add_reaction
|
|
|
|
!================================================================================
|
|
|
|
subroutine dealias_rhs(n)
|
|
|
|
use m_io
|
|
use m_parameters
|
|
use m_work
|
|
use x_fftw
|
|
|
|
implicit none
|
|
|
|
integer :: i, j, k, n
|
|
real*8 :: wnum2, akmax
|
|
|
|
akmax = real(kmax,8)
|
|
|
|
do k = 1,nz
|
|
do j = 1,ny
|
|
do i = 1,nx+1,2
|
|
|
|
if ( abs(akx(i)).gt.akmax .or. &
|
|
abs(aky(k)).gt.akmax .or. &
|
|
abs(akz(j)).gt.akmax ) then
|
|
|
|
wrk(i ,j,k,n) = zip
|
|
wrk(i+1,j,k,n) = zip
|
|
end if
|
|
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
return
|
|
|
|
end subroutine dealias_rhs
|
|
|
|
!================================================================================
|
|
!================================================================================
|
|
!================================================================================
|
|
!================================================================================
|
|
subroutine test_rhs_scalars
|
|
|
|
use m_openmpi
|
|
use m_io
|
|
use m_parameters
|
|
use m_fields
|
|
use m_work
|
|
use x_fftw
|
|
|
|
implicit none
|
|
|
|
integer :: i,j,k, n
|
|
real*8 :: a,b,c, x,y,z
|
|
if (task.eq.'hydro') then
|
|
|
|
a = 1.d0
|
|
b = 5.d0
|
|
c = 17.d0
|
|
|
|
do k = 1,nz
|
|
do j = 1,ny
|
|
do i = 1,nx
|
|
|
|
x = dx*real(i-1)
|
|
y = dx*real(j-1)
|
|
z = dx*real(myid*nz + k-1)
|
|
|
|
wrk(i,j,k,1) = sin(a * x)
|
|
wrk(i,j,k,2) = sin(b * y)
|
|
wrk(i,j,k,3) = sin(c * z)
|
|
wrk(i,j,k,4) = cos(a * x)
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
|
|
do n = 1,4
|
|
call xFFT3d(1,n)
|
|
fields(:,:,:,n) = wrk(:,:,:,n)
|
|
end do
|
|
|
|
|
|
nu = .5d0
|
|
|
|
call rhs_scalars
|
|
|
|
print *,'got rhs'
|
|
|
|
call xFFT3d(-1,4)
|
|
|
|
|
|
do k = 1,nz
|
|
do j = 1,ny
|
|
do i = 1,nx
|
|
|
|
x = dx*real(i-1)
|
|
y = dx*real(j-1)
|
|
z = dx*real(myid*nz + k-1)
|
|
|
|
|
|
! checking
|
|
wrk(i,j,k,0) = -a*cos(2.*a*x) - cos(a*x)*(b*cos(b*y) + c*cos(c*z) + nu*a**2)
|
|
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
!!$ tmp4(:,:,:) = wrk(1:nx,:,:,4)
|
|
!!$ fname = 'r1.arr'
|
|
!!$ call write_tmp4
|
|
!!$
|
|
!!$ tmp4(:,:,:) = wrk(1:nx,:,:,0)
|
|
!!$ fname = 'r0.arr'
|
|
!!$ call write_tmp4
|
|
|
|
|
|
|
|
|
|
wrk(:,:,:,0) = abs(wrk(1:nx,:,:,0) - wrk(1:nx,:,:,4))
|
|
|
|
print *,'Maximum error is ',maxval(wrk(1:nx,:,:,0))
|
|
|
|
!!$ tmp4(:,:,:) = wrk(1:nx,:,:,3) - wrk(1:nx,:,:,6)
|
|
!!$ fname = 'e3.arr'
|
|
!!$ call write_tmp4
|
|
!!$
|
|
|
|
|
|
end if
|
|
return
|
|
end subroutine test_rhs_scalars
|