gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Allow for gradual introduction of the forcing term.

Browse Source 
- the new parameter tbfor determines the time when the forcing starts
   to be introduced gradually; the parameter tefor determines at what
   time the forcing operates with the full power; the transition between
   tbfor and tefor is described by sinus function;
This commit is contained in:
Grzegorz Kowal 2011-06-10 15:46:40 -03:00
parent 3af808978d
commit 1abbafed99
4 changed files with 110 additions and 44 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
src/config.F90
View File

@ -133,6 +133,8 @@ module config
real , save :: ku = 3.0d+0 ! the maximum cut-off wave number
real , save :: kc = 0.1d+0 ! the spectrum width
integer , save :: kd = 1 ! the wave number increment
real , save :: tbfor = 0.0d+0 ! the moment when the forcing is initiated
real , save :: tefor = 0.0d+0 ! the moment when the forcing is reached full power
#endif /* FORCE */
#ifdef VISCOSITY
! viscous terms parameters
@ -389,6 +391,10 @@ module config
read(value,*) kc
case("kd")
read(value,*) kd
case("tbfor")
read(value,*) tbfor
case("tefor")
read(value,*) tefor
#endif /* FORCE */
#ifdef VISCOSITY
case("visc")

1
src/constants.F90
View File

@ -31,6 +31,7 @@ module constants
! number Pi and its multiplications
!
real, parameter :: hpi = 1.5707963267948965579989817342721d0
real, parameter :: pi = 3.1415926535897931159979634685442d0
real, parameter :: dpi = 6.2831853071795862319959269370884d0
real, parameter :: qpi = 12.5663706143591724639918538741770d0

122
src/evolution.F90
View File

@ -54,6 +54,7 @@ module evolution
use mesh , only : update_mesh
use timer , only : start_timer, stop_timer
#ifdef FORCE
use config , only : tbfor
use forcing , only : fourier_transform, evolve_forcing
use variables , only : idn, imz
#endif /* FORCE */
@ -72,23 +73,29 @@ module evolution
call start_timer(2)
#ifdef FORCE
! perform forcing evolution only when t >= tbfor
!
if (t .ge. tbfor) then
! perform the Fourier transform of the velocity field
!
pblock => list_data
do while (associated(pblock))
pblock => list_data
do while (associated(pblock))
if (pblock%meta%leaf) &
call fourier_transform(pblock%meta%level &
, pblock%meta%xmin, pblock%meta%ymin, pblock%meta%zmin &
, pblock%u(idn:imz,:,:,:))
if (pblock%meta%leaf) &
call fourier_transform(pblock%meta%level &
, pblock%meta%xmin, pblock%meta%ymin &
, pblock%meta%zmin, pblock%u(idn:imz,:,:,:))
pblock => pblock%next ! assign pointer to the next block
pblock => pblock%next ! assign pointer to the next block
end do
end do
! evolve the forcing source terms by the time interval dt
!
call evolve_forcing(dt)
call evolve_forcing(t, dt)
end if ! t >= tbfor
#endif /* FORCE */
! iterate over all data blocks and perform one step of time evolution
@ -333,6 +340,7 @@ module evolution
use variables, only : iph
#endif /* MHD & GLM */
#ifdef FORCE
use config , only : tbfor
use forcing , only : real_forcing
use variables, only : inx, iny, inz
use variables, only : idn, imx, imy, imz
@ -383,39 +391,45 @@ module evolution
!
pblock%u(iph,:,:,:) = decay * pblock%u(iph,:,:,:)
#endif /* MHD & GLM */
#ifdef FORCE
! add forcing term only if t >= tbfor
!
if (t .ge. tbfor) then
! obtain the forcing terms in real space
!
call real_forcing(pblock%meta%level, pblock%meta%xmin, pblock%meta%ymin &
, pblock%meta%zmin, f(:,:,:,:))
call real_forcing(pblock%meta%level, pblock%meta%xmin, pblock%meta%ymin &
, pblock%meta%zmin, f(:,:,:,:))
#ifdef ADI
! calculate kinetic energy before adding the forcing term
!
ek(:,:,:) = 0.5d0 * (pblock%u(imx,:,:,:)**2 + pblock%u(imy,:,:,:)**2 &
ek(:,:,:) = 0.5d0 * (pblock%u(imx,:,:,:)**2 + pblock%u(imy,:,:,:)**2 &
+ pblock%u(imz,:,:,:)**2) / pblock%u(idn,:,:,:)
#endif /* ADI */
! update momenta due to the forcing terms
!
pblock%u(imx,:,:,:) = pblock%u(imx,:,:,:) &
pblock%u(imx,:,:,:) = pblock%u(imx,:,:,:) &
+ pblock%u(idn,:,:,:) * f(inx,:,:,:)
pblock%u(imy,:,:,:) = pblock%u(imy,:,:,:) &
pblock%u(imy,:,:,:) = pblock%u(imy,:,:,:) &
+ pblock%u(idn,:,:,:) * f(iny,:,:,:)
pblock%u(imz,:,:,:) = pblock%u(imz,:,:,:) &
pblock%u(imz,:,:,:) = pblock%u(imz,:,:,:) &
+ pblock%u(idn,:,:,:) * f(inz,:,:,:)
#ifdef ADI
! calculate kinetic energy after adding the forcing term
!
dek(:,:,:) = 0.5d0 * (pblock%u(imx,:,:,:)**2 + pblock%u(imy,:,:,:)**2 &
dek(:,:,:) = 0.5d0 * (pblock%u(imx,:,:,:)**2 + pblock%u(imy,:,:,:)**2 &
+ pblock%u(imz,:,:,:)**2) / pblock%u(idn,:,:,:) - ek(:,:,:)
! update total energy with the injected one
!
pblock%u(ien,:,:,:) = pblock%u(ien,:,:,:) + dek(:,:,:)
pblock%u(ien,:,:,:) = pblock%u(ien,:,:,:) + dek(:,:,:)
#endif /* ADI */
end if ! t >= tbfor
#endif /* FORCE */
#ifdef SHAPE
@ -898,6 +912,7 @@ module evolution
use blocks , only : block_data
use config , only : im, jm, km
#ifdef FORCE
use config , only : tbfor
use forcing , only : real_forcing
#endif /* FORCE */
use coords , only : adxi, adyi, adzi
@ -975,35 +990,45 @@ module evolution
pblock%u(iph,:,:,:) = decay * pblock%u(iph,:,:,:)
#endif /* GLM */
#endif /* MHD */
#ifdef FORCE
! add the forcing term only if t >= tbfor
!
if (t .ge. tbfor) then
! obtain the forcing terms in real space
!
call real_forcing(pblock%meta%level, pblock%meta%xmin, pblock%meta%ymin &
call real_forcing(pblock%meta%level, pblock%meta%xmin, pblock%meta%ymin &
, pblock%meta%zmin, f(:,:,:,:))
#ifdef ADI
! calculate kinetic energy before adding the forcing term
!
ek(:,:,:) = 0.5d0 * (pblock%u(imx,:,:,:)**2 + pblock%u(imy,:,:,:)**2 &
ek(:,:,:) = 0.5d0 * (pblock%u(imx,:,:,:)**2 + pblock%u(imy,:,:,:)**2 &
+ pblock%u(imz,:,:,:)**2) / pblock%u(idn,:,:,:)
#endif /* ADI */
! update momenta due to the forcing terms
!
pblock%u(imx,:,:,:) = pblock%u(imx,:,:,:) + pblock%u(idn,:,:,:) * f(1,:,:,:)
pblock%u(imy,:,:,:) = pblock%u(imy,:,:,:) + pblock%u(idn,:,:,:) * f(2,:,:,:)
pblock%u(imz,:,:,:) = pblock%u(imz,:,:,:) + pblock%u(idn,:,:,:) * f(3,:,:,:)
pblock%u(imx,:,:,:) = pblock%u(imx,:,:,:) &
+ pblock%u(idn,:,:,:) * f(1,:,:,:)
pblock%u(imy,:,:,:) = pblock%u(imy,:,:,:) &
+ pblock%u(idn,:,:,:) * f(2,:,:,:)
pblock%u(imz,:,:,:) = pblock%u(imz,:,:,:) &
+ pblock%u(idn,:,:,:) * f(3,:,:,:)
#ifdef ADI
! calculate kinetic energy after adding the forcing term
!
dek(:,:,:) = 0.5d0 * (pblock%u(imx,:,:,:)**2 + pblock%u(imy,:,:,:)**2 &
dek(:,:,:) = 0.5d0 * (pblock%u(imx,:,:,:)**2 + pblock%u(imy,:,:,:)**2 &
+ pblock%u(imz,:,:,:)**2) / pblock%u(idn,:,:,:) - ek(:,:,:)
! update total energy with the injected one
!
pblock%u(ien,:,:,:) = pblock%u(ien,:,:,:) + dek(:,:,:)
pblock%u(ien,:,:,:) = pblock%u(ien,:,:,:) + dek(:,:,:)
#endif /* ADI */
end if ! t >= tbfor
#endif /* FORCE */
#ifdef SHAPE
@ -1030,6 +1055,7 @@ module evolution
use blocks , only : block_data
use config , only : im, jm, km
#ifdef FORCE
use config , only : tbfor
use forcing , only : real_forcing
#endif /* FORCE */
use coords , only : adxi, adyi, adzi
@ -1131,34 +1157,44 @@ module evolution
#endif /* GLM */
#endif /* MHD */
#ifdef FORCE
! add the forcing term only if t >= tbfor
!
if (t .ge. tbfor) then
! obtain the forcing terms in real space
!
call real_forcing(pblock%meta%level, pblock%meta%xmin, pblock%meta%ymin &
call real_forcing(pblock%meta%level, pblock%meta%xmin, pblock%meta%ymin &
, pblock%meta%zmin, f(:,:,:,:))
#ifdef ADI
! calculate kinetic energy before adding the forcing term
!
ek(:,:,:) = 0.5d0 * (pblock%u(imx,:,:,:)**2 + pblock%u(imy,:,:,:)**2 &
ek(:,:,:) = 0.5d0 * (pblock%u(imx,:,:,:)**2 + pblock%u(imy,:,:,:)**2 &
+ pblock%u(imz,:,:,:)**2) / pblock%u(idn,:,:,:)
#endif /* ADI */
! update momenta due to the forcing terms
!
pblock%u(imx,:,:,:) = pblock%u(imx,:,:,:) + pblock%u(idn,:,:,:) * f(1,:,:,:)
pblock%u(imy,:,:,:) = pblock%u(imy,:,:,:) + pblock%u(idn,:,:,:) * f(2,:,:,:)
pblock%u(imz,:,:,:) = pblock%u(imz,:,:,:) + pblock%u(idn,:,:,:) * f(3,:,:,:)
pblock%u(imx,:,:,:) = pblock%u(imx,:,:,:) &
+ pblock%u(idn,:,:,:) * f(1,:,:,:)
pblock%u(imy,:,:,:) = pblock%u(imy,:,:,:) &
+ pblock%u(idn,:,:,:) * f(2,:,:,:)
pblock%u(imz,:,:,:) = pblock%u(imz,:,:,:) &
+ pblock%u(idn,:,:,:) * f(3,:,:,:)
#ifdef ADI
! calculate kinetic energy after adding the forcing term
!
dek(:,:,:) = 0.5d0 * (pblock%u(imx,:,:,:)**2 + pblock%u(imy,:,:,:)**2 &
dek(:,:,:) = 0.5d0 * (pblock%u(imx,:,:,:)**2 + pblock%u(imy,:,:,:)**2 &
+ pblock%u(imz,:,:,:)**2) / pblock%u(idn,:,:,:) - ek(:,:,:)
! update total energy with the injected one
!
pblock%u(ien,:,:,:) = pblock%u(ien,:,:,:) + dek(:,:,:)
pblock%u(ien,:,:,:) = pblock%u(ien,:,:,:) + dek(:,:,:)
end if ! t >= tbfor
#endif /* ADI */
#endif /* FORCE */
@ -1187,6 +1223,7 @@ module evolution
use blocks , only : block_data
use config , only : im, jm, km
#ifdef FORCE
use config , only : tbfor
use forcing , only : real_forcing
#endif /* FORCE */
use coords , only : adxi, adyi, adzi
@ -1319,33 +1356,42 @@ module evolution
#endif /* MHD */
#ifdef FORCE
! add the forcing term only if t >= tbfor
!
if (t .ge. tbfor) then
! obtain the forcing terms in real space
!
call real_forcing(pblock%meta%level, pblock%meta%xmin, pblock%meta%ymin &
call real_forcing(pblock%meta%level, pblock%meta%xmin, pblock%meta%ymin &
, pblock%meta%zmin, f(:,:,:,:))
#ifdef ADI
! calculate kinetic energy before adding the forcing term
!
ek(:,:,:) = 0.5d0 * (pblock%u(imx,:,:,:)**2 + pblock%u(imy,:,:,:)**2 &
ek(:,:,:) = 0.5d0 * (pblock%u(imx,:,:,:)**2 + pblock%u(imy,:,:,:)**2 &
+ pblock%u(imz,:,:,:)**2) / pblock%u(idn,:,:,:)
#endif /* ADI */
! update momenta due to the forcing terms
!
pblock%u(imx,:,:,:) = pblock%u(imx,:,:,:) + pblock%u(idn,:,:,:) * f(1,:,:,:)
pblock%u(imy,:,:,:) = pblock%u(imy,:,:,:) + pblock%u(idn,:,:,:) * f(2,:,:,:)
pblock%u(imz,:,:,:) = pblock%u(imz,:,:,:) + pblock%u(idn,:,:,:) * f(3,:,:,:)
pblock%u(imx,:,:,:) = pblock%u(imx,:,:,:) &
+ pblock%u(idn,:,:,:) * f(1,:,:,:)
pblock%u(imy,:,:,:) = pblock%u(imy,:,:,:) &
+ pblock%u(idn,:,:,:) * f(2,:,:,:)
pblock%u(imz,:,:,:) = pblock%u(imz,:,:,:) &
+ pblock%u(idn,:,:,:) * f(3,:,:,:)
#ifdef ADI
! calculate kinetic energy after adding the forcing term
!
dek(:,:,:) = 0.5d0 * (pblock%u(imx,:,:,:)**2 + pblock%u(imy,:,:,:)**2 &
dek(:,:,:) = 0.5d0 * (pblock%u(imx,:,:,:)**2 + pblock%u(imy,:,:,:)**2 &
+ pblock%u(imz,:,:,:)**2) / pblock%u(idn,:,:,:) - ek(:,:,:)
! update total energy with the injected one
!
pblock%u(ien,:,:,:) = pblock%u(ien,:,:,:) + dek(:,:,:)
pblock%u(ien,:,:,:) = pblock%u(ien,:,:,:) + dek(:,:,:)
end if ! t >= tbfor
#endif /* ADI */
#endif /* FORCE */

25
src/forcing.F90
View File

@ -372,10 +372,10 @@ module forcing
!
!===============================================================================
!
subroutine evolve_forcing(dt)
subroutine evolve_forcing(t, dt)
use config , only : fpow, fdt
use constants, only : dpi
use config , only : fpow, fdt, tbfor, tefor
use constants, only : hpi, dpi
#ifdef MPI
use mpitools , only : mallreducesumc
#endif /* MPI */
@ -386,13 +386,14 @@ module forcing
! input arguments
!
real, intent(in) :: dt
real, intent(in) :: t, dt
! local variables
!
integer :: l, n, ni
complex :: aran, bran, xi1, xi2
real :: phi, psi, th1, th2, div, tanth1, fpw, fac
real :: ts, fc, ft
!-------------------------------------------------------------------------------
!
@ -405,8 +406,20 @@ module forcing
! reduce velocity fourier components from all processors
!
call mallreducesumc(nf, 3, vtab(:,:))
#endif /* MPI */
! calculate the amplitude increase factor
!
if (t .gt. tefor) then
ft = 1.0d0
else if ((t + dt) .lt. tbfor) then
ft = 0.0d0
else
ts = t - tbfor + 0.5d0 * dt
fc = hpi / max(epsilon(t), (tefor - tbfor))
ft = sin(fc * ts)
end if
! calculate the number of forcing integration iteration for the current timestep
!
ni = int(dt / fdt)
@ -471,7 +484,7 @@ module forcing
! normalize amplitudes to the correct input power
!
fpw = abs(sum(ftab(:,:) * conjg(ftab(:,:))))
fac = sqrt(4.0d0 * dt * fpow / fpw)
fac = ft * sqrt(4.0d0 * dt * fpow / fpw)
ftab(:,:) = fac * ftab(:,:)
! calculate the force-force correlation in Fourier space