Merge branch 'master' into reconnection
This commit is contained in:
commit
198fda648a
@ -31,10 +31,22 @@
|
||||
!
|
||||
module evolution
|
||||
|
||||
#ifdef PROFILE
|
||||
! import external subroutines
|
||||
!
|
||||
use timers, only : set_timer, start_timer, stop_timer
|
||||
#endif /* PROFILE */
|
||||
|
||||
! module variables are not implicit by default
|
||||
!
|
||||
implicit none
|
||||
|
||||
#ifdef PROFILE
|
||||
! timer indices
|
||||
!
|
||||
integer , save :: imi, ima, imt, imu, imf, imn, imv
|
||||
#endif /* PROFILE */
|
||||
|
||||
! pointer to the temporal integration subroutine
|
||||
!
|
||||
procedure(evolve_euler), pointer, save :: evolve => null()
|
||||
@ -116,6 +128,22 @@ module evolution
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
#ifdef PROFILE
|
||||
! set timer descriptions
|
||||
!
|
||||
call set_timer('evolution:: initialization' , imi)
|
||||
call set_timer('evolution:: solution advance', ima)
|
||||
call set_timer('evolution:: new time step' , imt)
|
||||
call set_timer('evolution:: solution update' , imu)
|
||||
call set_timer('evolution:: flux update' , imf)
|
||||
call set_timer('evolution:: increment update', imn)
|
||||
call set_timer('evolution:: variable update' , imv)
|
||||
|
||||
! start accounting time for module initialization/finalization
|
||||
!
|
||||
call start_timer(imi)
|
||||
#endif /* PROFILE */
|
||||
|
||||
! get the integration method and the value of the CFL coefficient
|
||||
!
|
||||
call get_parameter_string ("time_advance", integration)
|
||||
@ -185,6 +213,12 @@ module evolution
|
||||
|
||||
end if
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for module initialization/finalization
|
||||
!
|
||||
call stop_timer(imi)
|
||||
#endif /* PROFILE */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine initialize_evolution
|
||||
@ -213,9 +247,23 @@ module evolution
|
||||
integer, intent(inout) :: iret
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
#ifdef PROFILE
|
||||
! start accounting time for module initialization/finalization
|
||||
!
|
||||
call start_timer(imi)
|
||||
#endif /* PROFILE */
|
||||
|
||||
! nullify pointer to integration subroutine
|
||||
!
|
||||
nullify(evolve)
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for module initialization/finalization
|
||||
!
|
||||
call stop_timer(imi)
|
||||
#endif /* PROFILE */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine finalize_evolution
|
||||
@ -256,6 +304,12 @@ module evolution
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
#ifdef PROFILE
|
||||
! start accounting time for solution advance
|
||||
!
|
||||
call start_timer(ima)
|
||||
#endif /* PROFILE */
|
||||
|
||||
! find new time step
|
||||
!
|
||||
call new_time_step(dtnext)
|
||||
@ -296,6 +350,12 @@ module evolution
|
||||
|
||||
end if ! toplev > 1
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for solution advance
|
||||
!
|
||||
call stop_timer(ima)
|
||||
#endif /* PROFILE */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine advance
|
||||
@ -334,18 +394,17 @@ module evolution
|
||||
!
|
||||
implicit none
|
||||
|
||||
! input variables
|
||||
! subroutine arguments
|
||||
!
|
||||
real(kind=8), intent(in) :: dtnext
|
||||
real(kind=8), intent(in) :: dtnext
|
||||
|
||||
! local pointers
|
||||
!
|
||||
type(block_data), pointer :: pblock
|
||||
type(block_data), pointer :: pdata
|
||||
|
||||
! local variables
|
||||
!
|
||||
integer :: iret
|
||||
integer(kind=4) :: lev
|
||||
integer :: iret, mlev
|
||||
real(kind=8) :: cm, dx_min
|
||||
|
||||
! local parameters
|
||||
@ -354,53 +413,62 @@ module evolution
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
#ifdef PROFILE
|
||||
! start accounting time for new time step estimation
|
||||
!
|
||||
call start_timer(imt)
|
||||
#endif /* PROFILE */
|
||||
|
||||
! reset the maximum speed, and the highest level
|
||||
!
|
||||
cmax = eps
|
||||
lev = 1
|
||||
mlev = 1
|
||||
|
||||
! assign pdata with the first block on the data block list
|
||||
!
|
||||
pdata => list_data
|
||||
|
||||
! iterate over all data blocks in order to find the maximum speed among them
|
||||
! and the highest level which is required to obtain the spatial step
|
||||
! and the highest level which is required to obtain the minimum spacial step
|
||||
!
|
||||
pblock => list_data
|
||||
do while (associated(pblock))
|
||||
do while (associated(pdata))
|
||||
|
||||
! find the maximum level occupied by blocks (can be smaller than toplev)
|
||||
! update the maximum level
|
||||
!
|
||||
lev = max(lev, pblock%meta%level)
|
||||
mlev = max(mlev, pdata%meta%level)
|
||||
|
||||
! obtain the maximum speed for the current block
|
||||
! get the maximum speed for the current block
|
||||
!
|
||||
cm = maxspeed(pblock%q(:,:,:,:))
|
||||
cm = maxspeed(pdata%q(:,:,:,:))
|
||||
|
||||
! compare global and local maximum speeds
|
||||
! update the maximum speed
|
||||
!
|
||||
cmax = max(cmax, cm)
|
||||
|
||||
! assiociate the pointer with the next block
|
||||
! assign pdata to the next block
|
||||
!
|
||||
pblock => pblock%next
|
||||
pdata => pdata%next
|
||||
|
||||
end do
|
||||
end do ! over data blocks
|
||||
|
||||
#ifdef MPI
|
||||
! find maximum speed in the system from all processors
|
||||
! reduce maximum speed and level over all processors
|
||||
!
|
||||
call reduce_maximum_real (cmax, iret)
|
||||
call reduce_maximum_integer(lev , iret)
|
||||
call reduce_maximum_integer(mlev, iret)
|
||||
#endif /* MPI */
|
||||
|
||||
! calculate squared cmax
|
||||
! calculate the squared cmax
|
||||
!
|
||||
cmax2 = cmax * cmax
|
||||
|
||||
! find the smallest spatial step
|
||||
! find the smallest spacial step
|
||||
!
|
||||
#if NDIMS == 2
|
||||
dx_min = min(adx(lev), ady(lev))
|
||||
dx_min = min(adx(mlev), ady(mlev))
|
||||
#endif /* NDIMS == 2 */
|
||||
#if NDIMS == 3
|
||||
dx_min = min(adx(lev), ady(lev), adz(lev))
|
||||
dx_min = min(adx(mlev), ady(mlev), adz(mlev))
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
! calculate the new time step
|
||||
@ -408,13 +476,15 @@ module evolution
|
||||
dtn = cfl * dx_min / max(cmax &
|
||||
+ 2.0d+00 * max(viscosity, resistivity) / dx_min, eps)
|
||||
|
||||
! substitute the new time step
|
||||
!
|
||||
dt = dtn
|
||||
|
||||
! round the time
|
||||
!
|
||||
if (dtnext > 0.0d+00) dt = min(dt, dtnext)
|
||||
if (dtnext > 0.0d+00) dt = min(dtn, dtnext)
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for new time step estimation
|
||||
!
|
||||
call stop_timer(imt)
|
||||
#endif /* PROFILE */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
@ -434,6 +504,12 @@ module evolution
|
||||
! Subroutine advances the solution by one time step using the 1st order
|
||||
! Euler integration method.
|
||||
!
|
||||
! References:
|
||||
!
|
||||
! [1] Press, W. H, Teukolsky, S. A., Vetterling, W. T., Flannery, B. P.,
|
||||
! "Numerical Recipes in Fortran",
|
||||
! Cambridge University Press, Cambridge, 1992
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine evolve_euler()
|
||||
@ -455,7 +531,7 @@ module evolution
|
||||
|
||||
! local pointers
|
||||
!
|
||||
type(block_data), pointer :: pblock
|
||||
type(block_data), pointer :: pdata
|
||||
|
||||
! local arrays
|
||||
!
|
||||
@ -463,40 +539,51 @@ module evolution
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
! update fluxes for the first step of the RK2 integration
|
||||
#ifdef PROFILE
|
||||
! start accounting time for one step update
|
||||
!
|
||||
call start_timer(imu)
|
||||
#endif /* PROFILE */
|
||||
|
||||
! update fluxes
|
||||
!
|
||||
call update_fluxes()
|
||||
|
||||
! update the solution using numerical fluxes stored in the data blocks
|
||||
! assign pdata with the first block on the data block list
|
||||
!
|
||||
pblock => list_data
|
||||
do while (associated(pblock))
|
||||
pdata => list_data
|
||||
|
||||
! calculate variable increment for the current block
|
||||
! iterate over all data blocks
|
||||
!
|
||||
call update_increment(pblock, du(:,:,:,:))
|
||||
do while (associated(pdata))
|
||||
|
||||
! add source terms
|
||||
! calculate the variable increment
|
||||
!
|
||||
call update_sources(pblock, du(:,:,:,:))
|
||||
call update_increment(pdata, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! update the solution for the fluid variables
|
||||
! add the source terms
|
||||
!
|
||||
pblock%u0(1:nv,:,:,:) = pblock%u0(1:nv,:,:,:) + dt * du(1:nv,:,:,:)
|
||||
call update_sources(pdata, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! update the solution
|
||||
!
|
||||
pdata%u0(1:nv,1:im,1:jm,1:km) = pdata%u0(1:nv,1:im,1:jm,1:km) &
|
||||
+ dt * du(1:nv,1:im,1:jm,1:km)
|
||||
|
||||
! update the conservative variable pointer
|
||||
!
|
||||
pblock%u => pblock%u0
|
||||
pdata%u => pdata%u0
|
||||
|
||||
! update ψ by its source term
|
||||
! update ψ with its source term
|
||||
!
|
||||
if (ibp > 0) pblock%u(ibp,:,:,:) = decay * pblock%u(ibp,:,:,:)
|
||||
if (ibp > 0) pdata%u(ibp,1:im,1:jm,1:km) = &
|
||||
decay * pdata%u(ibp,1:im,1:jm,1:km)
|
||||
|
||||
! assign pointer to the next block
|
||||
! assign pdata to the next block
|
||||
!
|
||||
pblock => pblock%next
|
||||
pdata => pdata%next
|
||||
|
||||
end do
|
||||
end do ! over data blocks
|
||||
|
||||
! update primitive variables
|
||||
!
|
||||
@ -506,6 +593,12 @@ module evolution
|
||||
!
|
||||
call boundary_variables()
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for one step update
|
||||
!
|
||||
call stop_timer(imu)
|
||||
#endif /* PROFILE */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine evolve_euler
|
||||
@ -545,7 +638,7 @@ module evolution
|
||||
|
||||
! local pointers
|
||||
!
|
||||
type(block_data), pointer :: pblock
|
||||
type(block_data), pointer :: pdata
|
||||
|
||||
! local arrays
|
||||
!
|
||||
@ -553,36 +646,46 @@ module evolution
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
! update fluxes for the first step of the RK2 integration
|
||||
#ifdef PROFILE
|
||||
! start accounting time for one step update
|
||||
!
|
||||
call start_timer(imu)
|
||||
#endif /* PROFILE */
|
||||
|
||||
! update fluxes
|
||||
!
|
||||
call update_fluxes()
|
||||
|
||||
! update the solution using numerical fluxes stored in the data blocks
|
||||
! assign pdata with the first block on the data block list
|
||||
!
|
||||
pblock => list_data
|
||||
do while (associated(pblock))
|
||||
pdata => list_data
|
||||
|
||||
! calculate variable increment for the current block
|
||||
! iterate over all data blocks
|
||||
!
|
||||
call update_increment(pblock, du(:,:,:,:))
|
||||
do while (associated(pdata))
|
||||
|
||||
! add source terms
|
||||
! calculate the variable increment
|
||||
!
|
||||
call update_sources(pblock, du(:,:,:,:))
|
||||
call update_increment(pdata, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! update the solution for the fluid variables
|
||||
! add the source terms
|
||||
!
|
||||
pblock%u1(1:nv,:,:,:) = pblock%u0(1:nv,:,:,:) + dt * du(1:nv,:,:,:)
|
||||
call update_sources(pdata, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! update the intermediate solution
|
||||
!
|
||||
pdata%u1(1:nv,1:im,1:jm,1:km) = pdata%u0(1:nv,1:im,1:jm,1:km) &
|
||||
+ dt * du(1:nv,1:im,1:jm,1:km)
|
||||
|
||||
! update the conservative variable pointer
|
||||
!
|
||||
pblock%u => pblock%u1
|
||||
pdata%u => pdata%u1
|
||||
|
||||
! assign pointer to the next block
|
||||
! assign pdata to the next block
|
||||
!
|
||||
pblock => pblock%next
|
||||
pdata => pdata%next
|
||||
|
||||
end do
|
||||
end do ! over data blocks
|
||||
|
||||
! update primitive variables
|
||||
!
|
||||
@ -592,41 +695,46 @@ module evolution
|
||||
!
|
||||
call boundary_variables()
|
||||
|
||||
! update fluxes for the second step of the RK2 integration
|
||||
! update fluxes from the intermediate stage
|
||||
!
|
||||
call update_fluxes()
|
||||
|
||||
! update the solution using numerical fluxes stored in the data blocks
|
||||
! assign pdata with the first block on the data block list
|
||||
!
|
||||
pblock => list_data
|
||||
do while (associated(pblock))
|
||||
pdata => list_data
|
||||
|
||||
! calculate variable increment for the current block
|
||||
! iterate over all data blocks
|
||||
!
|
||||
call update_increment(pblock, du(:,:,:,:))
|
||||
do while (associated(pdata))
|
||||
|
||||
! add source terms
|
||||
! calculate the variable increment
|
||||
!
|
||||
call update_sources(pblock, du(:,:,:,:))
|
||||
call update_increment(pdata, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! update the solution for the fluid variables
|
||||
! add the source terms
|
||||
!
|
||||
pblock%u0(1:nv,:,:,:) = 0.5d+00 * (pblock%u0(1:nv,:,:,:) &
|
||||
+ pblock%u1(1:nv,:,:,:) + dt * du(1:nv,:,:,:))
|
||||
call update_sources(pdata, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! update the final solution
|
||||
!
|
||||
pdata%u0(1:nv,1:im,1:jm,1:km) = 0.5d+00 * (pdata%u0(1:nv,1:im,1:jm,1:km) &
|
||||
+ pdata%u1(1:nv,1:im,1:jm,1:km) &
|
||||
+ dt * du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! update the conservative variable pointer
|
||||
!
|
||||
pblock%u => pblock%u0
|
||||
pdata%u => pdata%u0
|
||||
|
||||
! update ψ by its source term
|
||||
! update ψ with its source term
|
||||
!
|
||||
if (ibp > 0) pblock%u(ibp,:,:,:) = decay * pblock%u(ibp,:,:,:)
|
||||
if (ibp > 0) pdata%u(ibp,1:im,1:jm,1:km) = &
|
||||
decay * pdata%u(ibp,1:im,1:jm,1:km)
|
||||
|
||||
! assign pointer to the next block
|
||||
! assign pdata to the next block
|
||||
!
|
||||
pblock => pblock%next
|
||||
pdata => pdata%next
|
||||
|
||||
end do
|
||||
end do ! over data blocks
|
||||
|
||||
! update primitive variables
|
||||
!
|
||||
@ -636,6 +744,12 @@ module evolution
|
||||
!
|
||||
call boundary_variables()
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for one step update
|
||||
!
|
||||
call stop_timer(imu)
|
||||
#endif /* PROFILE */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine evolve_rk2
|
||||
@ -696,6 +810,12 @@ module evolution
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
#ifdef PROFILE
|
||||
! start accounting time for one step update
|
||||
!
|
||||
call start_timer(imu)
|
||||
#endif /* PROFILE */
|
||||
|
||||
! prepare things which don't change later
|
||||
!
|
||||
if (first) then
|
||||
@ -819,7 +939,7 @@ module evolution
|
||||
!
|
||||
pdata%u => pdata%u0
|
||||
|
||||
! update ψ by its source term
|
||||
! update ψ with its source term
|
||||
!
|
||||
if (ibp > 0) pdata%u(ibp,1:im,1:jm,1:km) = &
|
||||
decay * pdata%u(ibp,1:im,1:jm,1:km)
|
||||
@ -838,6 +958,12 @@ module evolution
|
||||
!
|
||||
call boundary_variables()
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for one step update
|
||||
!
|
||||
call stop_timer(imu)
|
||||
#endif /* PROFILE */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine evolve_ssprk2
|
||||
@ -878,7 +1004,7 @@ module evolution
|
||||
|
||||
! local pointers
|
||||
!
|
||||
type(block_data), pointer :: pblock
|
||||
type(block_data), pointer :: pdata
|
||||
|
||||
! local variables
|
||||
!
|
||||
@ -895,42 +1021,52 @@ module evolution
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
#ifdef PROFILE
|
||||
! start accounting time for one step update
|
||||
!
|
||||
call start_timer(imu)
|
||||
#endif /* PROFILE */
|
||||
|
||||
!! 1st substep of integration
|
||||
!!
|
||||
! prepare fractional time step
|
||||
! prepare the fractional time step
|
||||
!
|
||||
ds = dt
|
||||
|
||||
! update fluxes for the first step of the RK2 integration
|
||||
! update fluxes
|
||||
!
|
||||
call update_fluxes()
|
||||
|
||||
! update the solution using numerical fluxes stored in the data blocks
|
||||
! assign pdata with the first block on the data block list
|
||||
!
|
||||
pblock => list_data
|
||||
do while (associated(pblock))
|
||||
pdata => list_data
|
||||
|
||||
! calculate variable increment for the current block
|
||||
! iterate over all data blocks
|
||||
!
|
||||
call update_increment(pblock, du(:,:,:,:))
|
||||
do while (associated(pdata))
|
||||
|
||||
! add source terms
|
||||
! calculate the variable increment
|
||||
!
|
||||
call update_sources(pblock, du(:,:,:,:))
|
||||
call update_increment(pdata, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! update the solution for the fluid variables
|
||||
! add the source terms
|
||||
!
|
||||
pblock%u1(1:nv,:,:,:) = pblock%u0(1:nv,:,:,:) + ds * du(1:nv,:,:,:)
|
||||
call update_sources(pdata, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! update the first intermediate solution
|
||||
!
|
||||
pdata%u1(1:nv,1:im,1:jm,1:km) = pdata%u0(1:nv,1:im,1:jm,1:km) &
|
||||
+ ds * du(1:nv,1:im,1:jm,1:km)
|
||||
|
||||
! update the conservative variable pointer
|
||||
!
|
||||
pblock%u => pblock%u1
|
||||
pdata%u => pdata%u1
|
||||
|
||||
! assign pointer to the next block
|
||||
! assign pdata to the next block
|
||||
!
|
||||
pblock => pblock%next
|
||||
pdata => pdata%next
|
||||
|
||||
end do
|
||||
end do ! over data blocks
|
||||
|
||||
! update primitive variables
|
||||
!
|
||||
@ -942,37 +1078,41 @@ module evolution
|
||||
|
||||
!! 2nd substep of integration
|
||||
!!
|
||||
! prepare fractional time step
|
||||
! prepare the fractional time step
|
||||
!
|
||||
ds = f22 * dt
|
||||
|
||||
! update fluxes for the first step of the RK2 integration
|
||||
! update fluxes
|
||||
!
|
||||
call update_fluxes()
|
||||
|
||||
! update the solution using numerical fluxes stored in the data blocks
|
||||
! assign pdata with the first block on the data block list
|
||||
!
|
||||
pblock => list_data
|
||||
do while (associated(pblock))
|
||||
pdata => list_data
|
||||
|
||||
! calculate variable increment for the current block
|
||||
! iterate over all data blocks
|
||||
!
|
||||
call update_increment(pblock, du(:,:,:,:))
|
||||
do while (associated(pdata))
|
||||
|
||||
! add source terms
|
||||
! calculate the variable increment
|
||||
!
|
||||
call update_sources(pblock, du(:,:,:,:))
|
||||
call update_increment(pdata, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! update the solution for the fluid variables
|
||||
! add the source terms
|
||||
!
|
||||
pblock%u1(1:nv,:,:,:) = f21 * pblock%u0(1:nv,:,:,:) &
|
||||
+ f22 * pblock%u1(1:nv,:,:,:) + ds * du(1:nv,:,:,:)
|
||||
call update_sources(pdata, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! assign pointer to the next block
|
||||
! update the second intermediate solution
|
||||
!
|
||||
pblock => pblock%next
|
||||
pdata%u1(1:nv,1:im,1:jm,1:km) = f21 * pdata%u0(1:nv,1:im,1:jm,1:km) &
|
||||
+ f22 * pdata%u1(1:nv,1:im,1:jm,1:km) &
|
||||
+ ds * du(1:nv,1:im,1:jm,1:km)
|
||||
|
||||
end do
|
||||
! assign pdata to the next block
|
||||
!
|
||||
pdata => pdata%next
|
||||
|
||||
end do ! over data blocks
|
||||
|
||||
! update primitive variables
|
||||
!
|
||||
@ -984,45 +1124,50 @@ module evolution
|
||||
|
||||
!! 3rd substep of integration
|
||||
!!
|
||||
! prepare fractional time step
|
||||
! prepare the fractional time step
|
||||
!
|
||||
ds = f32 * dt
|
||||
|
||||
! update fluxes for the second step of the RK2 integration
|
||||
! update fluxes
|
||||
!
|
||||
call update_fluxes()
|
||||
|
||||
! update the solution using numerical fluxes stored in the data blocks
|
||||
! assign pdata with the first block on the data block list
|
||||
!
|
||||
pblock => list_data
|
||||
do while (associated(pblock))
|
||||
pdata => list_data
|
||||
|
||||
! calculate variable increment for the current block
|
||||
! iterate over all data blocks
|
||||
!
|
||||
call update_increment(pblock, du(:,:,:,:))
|
||||
do while (associated(pdata))
|
||||
|
||||
! add source terms
|
||||
! calculate the variable increment
|
||||
!
|
||||
call update_sources(pblock, du(:,:,:,:))
|
||||
call update_increment(pdata, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! update the solution for the fluid variables
|
||||
! add the source terms
|
||||
!
|
||||
pblock%u0(1:nv,:,:,:) = f31 * pblock%u0(1:nv,:,:,:) &
|
||||
+ f32 * pblock%u1(1:nv,:,:,:) + ds * du(1:nv,:,:,:)
|
||||
call update_sources(pdata, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! update the final solution
|
||||
!
|
||||
pdata%u0(1:nv,1:im,1:jm,1:km) = f31 * pdata%u0(1:nv,1:im,1:jm,1:km) &
|
||||
+ f32 * pdata%u1(1:nv,1:im,1:jm,1:km) &
|
||||
+ ds * du(1:nv,1:im,1:jm,1:km)
|
||||
|
||||
! update the conservative variable pointer
|
||||
!
|
||||
pblock%u => pblock%u0
|
||||
pdata%u => pdata%u0
|
||||
|
||||
! update ψ by its source term
|
||||
! update ψ with its source term
|
||||
!
|
||||
if (ibp > 0) pblock%u(ibp,:,:,:) = decay * pblock%u(ibp,:,:,:)
|
||||
if (ibp > 0) pdata%u(ibp,1:im,1:jm,1:km) = &
|
||||
decay * pdata%u(ibp,1:im,1:jm,1:km)
|
||||
|
||||
! assign pointer to the next block
|
||||
! assign pdata to the next block
|
||||
!
|
||||
pblock => pblock%next
|
||||
pdata => pdata%next
|
||||
|
||||
end do
|
||||
end do ! over data blocks
|
||||
|
||||
! update primitive variables
|
||||
!
|
||||
@ -1032,6 +1177,12 @@ module evolution
|
||||
!
|
||||
call boundary_variables()
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for one step update
|
||||
!
|
||||
call stop_timer(imu)
|
||||
#endif /* PROFILE */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine evolve_rk3
|
||||
@ -1090,6 +1241,12 @@ module evolution
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
#ifdef PROFILE
|
||||
! start accounting time for one step update
|
||||
!
|
||||
call start_timer(imu)
|
||||
#endif /* PROFILE */
|
||||
|
||||
!= 1st step: U(1) = U(n) + 1/2 dt F[U(n)]
|
||||
!
|
||||
! calculate the fractional time step
|
||||
@ -1261,7 +1418,7 @@ module evolution
|
||||
!
|
||||
pdata%u => pdata%u0
|
||||
|
||||
! update ψ by its source term
|
||||
! update ψ with its source term
|
||||
!
|
||||
if (ibp > 0) pdata%u(ibp,1:im,1:jm,1:km) = &
|
||||
decay * pdata%u(ibp,1:im,1:jm,1:km)
|
||||
@ -1280,6 +1437,12 @@ module evolution
|
||||
!
|
||||
call boundary_variables()
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for one step update
|
||||
!
|
||||
call stop_timer(imu)
|
||||
#endif /* PROFILE */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine evolve_ssprk34
|
||||
@ -1346,6 +1509,12 @@ module evolution
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
#ifdef PROFILE
|
||||
! start accounting time for one step update
|
||||
!
|
||||
call start_timer(imu)
|
||||
#endif /* PROFILE */
|
||||
|
||||
!= 1st step: U(1) = U(n) + b1 dt F[U(n)]
|
||||
!
|
||||
! calculate the fractional time step
|
||||
@ -1564,7 +1733,7 @@ module evolution
|
||||
+ a55 * pdata%u0(1:nv,1:im,1:jm,1:km) &
|
||||
+ ds * du(1:nv,1:im,1:jm,1:km)
|
||||
|
||||
! update ψ by its source term
|
||||
! update ψ with its source term
|
||||
!
|
||||
if (ibp > 0) pdata%u(ibp,1:im,1:jm,1:km) = &
|
||||
decay * pdata%u(ibp,1:im,1:jm,1:km)
|
||||
@ -1583,6 +1752,12 @@ module evolution
|
||||
!
|
||||
call boundary_variables()
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for one step update
|
||||
!
|
||||
call stop_timer(imu)
|
||||
#endif /* PROFILE */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine evolve_ssprk35
|
||||
@ -1610,6 +1785,8 @@ module evolution
|
||||
!
|
||||
use blocks , only : block_data, list_data
|
||||
use coordinates , only : adx, ady, adz
|
||||
use coordinates , only : im, jm, km
|
||||
use equations , only : nv
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
@ -1617,7 +1794,7 @@ module evolution
|
||||
|
||||
! local pointers
|
||||
!
|
||||
type(block_data), pointer :: pblock
|
||||
type(block_data), pointer :: pdata
|
||||
|
||||
! local vectors
|
||||
!
|
||||
@ -1629,34 +1806,50 @@ module evolution
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
#ifdef PROFILE
|
||||
! start accounting time for fluxe update
|
||||
!
|
||||
call start_timer(imf)
|
||||
#endif /* PROFILE */
|
||||
|
||||
! assign pdata with the first block on the data block list
|
||||
!
|
||||
pdata => list_data
|
||||
|
||||
! iterate over all data blocks
|
||||
!
|
||||
pblock => list_data
|
||||
do while (associated(pblock))
|
||||
do while (associated(pdata))
|
||||
|
||||
! obtain dx, dy, and dz for the current block
|
||||
!
|
||||
dx(1) = adx(pblock%meta%level)
|
||||
dx(2) = ady(pblock%meta%level)
|
||||
dx(3) = adz(pblock%meta%level)
|
||||
dx(1) = adx(pdata%meta%level)
|
||||
dx(2) = ady(pdata%meta%level)
|
||||
dx(3) = adz(pdata%meta%level)
|
||||
|
||||
! update the flux for the current block
|
||||
! update fluxes for the current block
|
||||
!
|
||||
do n = 1, NDIMS
|
||||
call update_flux(n, dx(n), pblock%q(:,:,:,:), pblock%f(n,:,:,:,:))
|
||||
call update_flux(n, dx(n), pdata%q(1:nv,1:im,1:jm,1:km) &
|
||||
, pdata%f(n,1:nv,1:im,1:jm,1:km))
|
||||
end do
|
||||
|
||||
! assign pointer to the next block
|
||||
! assign pdata to the next block
|
||||
!
|
||||
pblock => pblock%next
|
||||
pdata => pdata%next
|
||||
|
||||
end do
|
||||
end do ! over data blocks
|
||||
|
||||
! correct the numerical fluxes of the blocks which have neighbours at higher
|
||||
! level
|
||||
! levels
|
||||
!
|
||||
call boundary_fluxes()
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for flux update
|
||||
!
|
||||
call stop_timer(imf)
|
||||
#endif /* PROFILE */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine update_fluxes
|
||||
@ -1701,6 +1894,12 @@ module evolution
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
#ifdef PROFILE
|
||||
! start accounting time for increment update
|
||||
!
|
||||
call start_timer(imn)
|
||||
#endif /* PROFILE */
|
||||
|
||||
! reset the increment array du
|
||||
!
|
||||
du(:,:,:,:) = 0.0d+00
|
||||
@ -1736,6 +1935,12 @@ module evolution
|
||||
end do
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for increment update
|
||||
!
|
||||
call stop_timer(imn)
|
||||
#endif /* PROFILE */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine update_increment
|
||||
@ -1774,6 +1979,12 @@ module evolution
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
#ifdef PROFILE
|
||||
! start accounting time for variable update
|
||||
!
|
||||
call start_timer(imv)
|
||||
#endif /* PROFILE */
|
||||
|
||||
! associate the pointer with the first block on the data block list
|
||||
!
|
||||
pdata => list_data
|
||||
@ -1800,6 +2011,12 @@ module evolution
|
||||
|
||||
end do
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for variable update
|
||||
!
|
||||
call stop_timer(imv)
|
||||
#endif /* PROFILE */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine update_variables
|
||||
|
||||
@ -148,7 +148,7 @@ mesh.o : mesh.F90 blocks.o coordinates.o domains.o equations.o \
|
||||
error.o interpolations.o mpitools.o problems.o refinement.o \
|
||||
timers.o
|
||||
mpitools.o : mpitools.F90 timers.o
|
||||
operators.o : operators.F90
|
||||
operators.o : operators.F90 timers.o
|
||||
parameters.o : parameters.F90 mpitools.o
|
||||
problems.o : problems.F90 blocks.o constants.o coordinates.o equations.o \
|
||||
error.o operators.o parameters.o random.o timers.o
|
||||
|
||||
32
src/mesh.F90
32
src/mesh.F90
@ -878,7 +878,7 @@ module mesh
|
||||
|
||||
! local buffer for data block exchange
|
||||
!
|
||||
real(kind=8) , dimension(2,nv,im,jm,km) :: rbuf
|
||||
real(kind=8) , dimension(3,nv,im,jm,km) :: rbuf
|
||||
#endif /* MPI */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
@ -931,8 +931,9 @@ module mesh
|
||||
|
||||
! copy data to buffer
|
||||
!
|
||||
rbuf(1,:,:,:,:) = pmeta%data%u(:,:,:,:)
|
||||
rbuf(2,:,:,:,:) = pmeta%data%q(:,:,:,:)
|
||||
rbuf(1,:,:,:,:) = pmeta%data%q (:,:,:,:)
|
||||
rbuf(2,:,:,:,:) = pmeta%data%u0(:,:,:,:)
|
||||
rbuf(3,:,:,:,:) = pmeta%data%u1(:,:,:,:)
|
||||
|
||||
! send data
|
||||
!
|
||||
@ -961,8 +962,9 @@ module mesh
|
||||
|
||||
! coppy the buffer to data block
|
||||
!
|
||||
pmeta%data%u(:,:,:,:) = rbuf(1,:,:,:,:)
|
||||
pmeta%data%q(:,:,:,:) = rbuf(2,:,:,:,:)
|
||||
pmeta%data%q (:,:,:,:) = rbuf(1,:,:,:,:)
|
||||
pmeta%data%u0(:,:,:,:) = rbuf(2,:,:,:,:)
|
||||
pmeta%data%u1(:,:,:,:) = rbuf(3,:,:,:,:)
|
||||
|
||||
end if ! nproc == n
|
||||
|
||||
@ -1031,11 +1033,11 @@ module mesh
|
||||
!
|
||||
! Arguments:
|
||||
!
|
||||
! pblock - the input meta block;
|
||||
! pmeta - the input meta block;
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine prolong_block(pblock)
|
||||
subroutine prolong_block(pmeta)
|
||||
|
||||
! import external procedures and variables
|
||||
!
|
||||
@ -1051,7 +1053,7 @@ module mesh
|
||||
|
||||
! input arguments
|
||||
!
|
||||
type(block_meta), pointer, intent(inout) :: pblock
|
||||
type(block_meta), pointer, intent(inout) :: pmeta
|
||||
|
||||
! local variables
|
||||
!
|
||||
@ -1083,7 +1085,7 @@ module mesh
|
||||
|
||||
! assign the pdata pointer
|
||||
!
|
||||
pdata => pblock%data
|
||||
pdata => pmeta%data
|
||||
|
||||
! prepare dimensions
|
||||
!
|
||||
@ -1175,7 +1177,7 @@ module mesh
|
||||
|
||||
! assign pointer to the current child
|
||||
!
|
||||
pchild => pblock%child(p)%ptr
|
||||
pchild => pmeta%child(p)%ptr
|
||||
|
||||
! obtain the position of child in the parent block
|
||||
!
|
||||
@ -1236,11 +1238,11 @@ module mesh
|
||||
!
|
||||
! Arguments:
|
||||
!
|
||||
! pblock - the input meta block;
|
||||
! pmeta - the input meta block;
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine restrict_block(pblock)
|
||||
subroutine restrict_block(pmeta)
|
||||
|
||||
! import external procedures and variables
|
||||
!
|
||||
@ -1256,7 +1258,7 @@ module mesh
|
||||
|
||||
! subroutine arguments
|
||||
!
|
||||
type(block_meta), pointer, intent(inout) :: pblock
|
||||
type(block_meta), pointer, intent(inout) :: pmeta
|
||||
|
||||
! local variables
|
||||
!
|
||||
@ -1282,7 +1284,7 @@ module mesh
|
||||
|
||||
! assign the parent data pointer
|
||||
!
|
||||
pparent => pblock%data
|
||||
pparent => pmeta%data
|
||||
|
||||
! iterate over all children
|
||||
!
|
||||
@ -1290,7 +1292,7 @@ module mesh
|
||||
|
||||
! assign a pointer to the current child
|
||||
!
|
||||
pchild => pblock%child(p)%ptr%data
|
||||
pchild => pmeta%child(p)%ptr%data
|
||||
|
||||
! obtain the child position in the parent block
|
||||
!
|
||||
|
||||
@ -30,10 +30,22 @@
|
||||
!
|
||||
module operators
|
||||
|
||||
#ifdef PROFILE
|
||||
! import external subroutines
|
||||
!
|
||||
use timers, only : set_timer, start_timer, stop_timer
|
||||
#endif /* PROFILE */
|
||||
|
||||
! module variables are not implicit by default
|
||||
!
|
||||
implicit none
|
||||
|
||||
#ifdef PROFILE
|
||||
! timer indices
|
||||
!
|
||||
integer , save :: imi, imd, img, imc, iml, im1, im2
|
||||
#endif /* PROFILE */
|
||||
|
||||
! by default everything is public
|
||||
!
|
||||
private
|
||||
@ -41,7 +53,7 @@ module operators
|
||||
! declare public subroutines
|
||||
!
|
||||
public :: initialize_operators, finalize_operators
|
||||
public :: divergence, curl, laplace
|
||||
public :: divergence, gradient, curl, laplace
|
||||
|
||||
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
||||
!
|
||||
@ -80,6 +92,27 @@ module operators
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
#ifdef PROFILE
|
||||
! set timer descriptions
|
||||
!
|
||||
call set_timer('operators:: initialization', imi)
|
||||
call set_timer('operators:: divergence' , imd)
|
||||
call set_timer('operators:: gradient' , img)
|
||||
call set_timer('operators:: curl' , imc)
|
||||
call set_timer('operators:: laplace' , iml)
|
||||
call set_timer('operators:: 1st derivative', im1)
|
||||
call set_timer('operators:: 2nd derivative', im2)
|
||||
|
||||
! start accounting time for the module initialization/finalization
|
||||
!
|
||||
call start_timer(imi)
|
||||
#endif /* PROFILE */
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for the module initialization/finalization
|
||||
!
|
||||
call stop_timer(imi)
|
||||
#endif /* PROFILE */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
@ -110,6 +143,17 @@ module operators
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
#ifdef PROFILE
|
||||
! start accounting time for the module initialization/finalization
|
||||
!
|
||||
call start_timer(imi)
|
||||
#endif /* PROFILE */
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for the module initialization/finalization
|
||||
!
|
||||
call stop_timer(imi)
|
||||
#endif /* PROFILE */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
@ -155,6 +199,12 @@ module operators
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
#ifdef PROFILE
|
||||
! start accounting time for the divergence operator calculation
|
||||
!
|
||||
call start_timer(imd)
|
||||
#endif /* PROFILE */
|
||||
|
||||
! allocate temporary array
|
||||
!
|
||||
allocate(w(size(u,2), size(u,3), size(u,4)))
|
||||
@ -167,7 +217,7 @@ module operators
|
||||
!
|
||||
do dir = 1, NDIMS
|
||||
|
||||
! calculate contribution from the Y derivative of the Y component
|
||||
! calculate derivative along the current direction
|
||||
!
|
||||
call derivative_1st(dir, dh(dir), u(dir,:,:,:), w(:,:,:))
|
||||
|
||||
@ -181,12 +231,83 @@ module operators
|
||||
!
|
||||
deallocate(w)
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for the divergence operator calculation
|
||||
!
|
||||
call stop_timer(imd)
|
||||
#endif /* PROFILE */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine divergence
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine GRADIENT:
|
||||
! -------------------
|
||||
!
|
||||
! Subroutine calculates the cell centered gradient of the input scalar field.
|
||||
!
|
||||
! grad(U) = ∇ U = [ ∂x U, ∂y U, ∂z U ]
|
||||
!
|
||||
! Arguments:
|
||||
!
|
||||
! dh - the spacial intervals in all direction;
|
||||
! u - the input scalar field;
|
||||
! v - the output gradient vector field;
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine gradient(dh, u, v)
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
implicit none
|
||||
|
||||
! input and output variables
|
||||
!
|
||||
real(kind=8), dimension(3) , intent(in) :: dh
|
||||
real(kind=8), dimension(:,:,:) , intent(in) :: u
|
||||
real(kind=8), dimension(:,:,:,:), intent(out) :: v
|
||||
|
||||
! local variables
|
||||
!
|
||||
integer :: dir
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
#ifdef PROFILE
|
||||
! start accounting time for the gradient operator calculation
|
||||
!
|
||||
call start_timer(img)
|
||||
#endif /* PROFILE */
|
||||
|
||||
! reset the output array
|
||||
!
|
||||
v(:,:,:,:) = 0.0d+00
|
||||
|
||||
! iterate over directions and calculate gradient components
|
||||
!
|
||||
do dir = 1, NDIMS
|
||||
|
||||
! calculate derivative along the current direction
|
||||
!
|
||||
call derivative_1st(dir, dh(dir), u(:,:,:), v(dir,:,:,:))
|
||||
|
||||
end do ! directions
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for the gradient operator calculation
|
||||
!
|
||||
call stop_timer(img)
|
||||
#endif /* PROFILE */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine gradient
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine CURL:
|
||||
! ---------------
|
||||
!
|
||||
@ -221,6 +342,12 @@ module operators
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
#ifdef PROFILE
|
||||
! start accounting time for the rotation operator calculation
|
||||
!
|
||||
call start_timer(imc)
|
||||
#endif /* PROFILE */
|
||||
|
||||
! allocate temporary array
|
||||
!
|
||||
allocate(w(size(u,2), size(u,3), size(u,4)))
|
||||
@ -295,6 +422,12 @@ module operators
|
||||
!
|
||||
deallocate(w)
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for the rotation operator calculation
|
||||
!
|
||||
call stop_timer(imc)
|
||||
#endif /* PROFILE */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine curl
|
||||
@ -334,6 +467,12 @@ module operators
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
#ifdef PROFILE
|
||||
! start accounting time for the laplace operator calculation
|
||||
!
|
||||
call start_timer(iml)
|
||||
#endif /* PROFILE */
|
||||
|
||||
! allocate temporary array
|
||||
!
|
||||
allocate(w(size(u,1), size(u,2), size(u,3)))
|
||||
@ -368,6 +507,12 @@ module operators
|
||||
!
|
||||
deallocate(w)
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for the laplace operator calculation
|
||||
!
|
||||
call stop_timer(iml)
|
||||
#endif /* PROFILE */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine laplace
|
||||
@ -420,6 +565,12 @@ module operators
|
||||
if (dir < 1 .or. dir > NDIMS .or. dh == 0.0d+00) return
|
||||
#endif /* DEBUG */
|
||||
|
||||
#ifdef PROFILE
|
||||
! start accounting time for the 1st derivative calculation
|
||||
!
|
||||
call start_timer(im1)
|
||||
#endif /* PROFILE */
|
||||
|
||||
! prepare index limits
|
||||
!
|
||||
m0 = size(u, dir)
|
||||
@ -458,6 +609,12 @@ module operators
|
||||
|
||||
end select
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for the 1st derivative calculation
|
||||
!
|
||||
call stop_timer(im1)
|
||||
#endif /* PROFILE */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine derivative_1st
|
||||
@ -504,6 +661,12 @@ module operators
|
||||
if (dir < 1 .or. dir > NDIMS .or. dh == 0.0d+00) return
|
||||
#endif /* DEBUG */
|
||||
|
||||
#ifdef PROFILE
|
||||
! start accounting time for the 2nd derivative calculation
|
||||
!
|
||||
call start_timer(im2)
|
||||
#endif /* PROFILE */
|
||||
|
||||
! prepare index limits
|
||||
!
|
||||
m0 = size(u, dir)
|
||||
@ -542,6 +705,12 @@ module operators
|
||||
|
||||
end select
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for the 2nd derivative calculation
|
||||
!
|
||||
call stop_timer(im2)
|
||||
#endif /* PROFILE */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine derivative_2nd
|
||||
|
||||
129
src/sources.F90
129
src/sources.F90
@ -45,6 +45,10 @@ module sources
|
||||
integer, save :: imi, imu
|
||||
#endif /* PROFILE */
|
||||
|
||||
! GLM-MHD source terms type (1 - EGLM, 2 - HEGLM)
|
||||
!
|
||||
integer , save :: glm_type = 0
|
||||
|
||||
! gravitational acceleration coefficient
|
||||
!
|
||||
real(kind=8), save :: gpoint = 0.0d+00
|
||||
@ -95,7 +99,7 @@ module sources
|
||||
|
||||
! include external procedures and variables
|
||||
!
|
||||
use parameters , only : get_parameter_real
|
||||
use parameters , only : get_parameter_string, get_parameter_real
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
@ -105,6 +109,10 @@ module sources
|
||||
!
|
||||
logical, intent(in) :: verbose
|
||||
integer, intent(inout) :: iret
|
||||
|
||||
! local variables
|
||||
!
|
||||
character(len=8) :: tglm = "none"
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
@ -119,6 +127,23 @@ module sources
|
||||
call start_timer(imi)
|
||||
#endif /* PROFILE */
|
||||
|
||||
! get the type of the GLM source terms
|
||||
!
|
||||
call get_parameter_string("glm_source_terms", tglm)
|
||||
|
||||
! set the glm_type variable to correct value
|
||||
!
|
||||
select case(trim(tglm))
|
||||
case("eglm", "EGLM")
|
||||
glm_type = 1
|
||||
|
||||
case("heglm", "HEGLM")
|
||||
glm_type = 2
|
||||
|
||||
case default
|
||||
glm_type = 0
|
||||
end select
|
||||
|
||||
! get acceleration coefficient
|
||||
!
|
||||
call get_parameter_real("gpoint" , gpoint )
|
||||
@ -135,6 +160,7 @@ module sources
|
||||
!
|
||||
if (verbose) then
|
||||
|
||||
write (*,"(4x,a,1x,a) ") "glm source terms =", trim(tglm)
|
||||
write (*,"(4x,a,1x,1e9.2)") "point mass constant =", gpoint
|
||||
write (*,"(4x,a,1x,1e9.2)") "viscosity =", viscosity
|
||||
write (*,"(4x,a,1x,1e9.2)") "resistivity =", resistivity
|
||||
@ -215,8 +241,8 @@ module sources
|
||||
use coordinates , only : ax, ay, az, adx, ady, adz, adxi, adyi, adzi
|
||||
use equations , only : nv, inx, iny, inz
|
||||
use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ien
|
||||
use equations , only : ibx, iby, ibz
|
||||
use operators , only : laplace, curl
|
||||
use equations , only : ibx, iby, ibz, ibp
|
||||
use operators , only : divergence, gradient, laplace, curl
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
@ -241,6 +267,7 @@ module sources
|
||||
real(kind=8), dimension(im) :: x
|
||||
real(kind=8), dimension(jm) :: y
|
||||
real(kind=8), dimension(km) :: z
|
||||
real(kind=8), dimension(im,jm,km) :: db
|
||||
real(kind=8), dimension(3,im,jm,km) :: jc
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
@ -440,9 +467,9 @@ module sources
|
||||
|
||||
end if ! viscosity is not zero
|
||||
|
||||
! proceed only if the resistivity coefficient is not zero
|
||||
!=== add magnetic field related source terms ===
|
||||
!
|
||||
if (resistivity > 0.0d+00 .and. ibx > 0) then
|
||||
if (ibx > 0) then
|
||||
|
||||
! prepare coordinate increments
|
||||
!
|
||||
@ -450,43 +477,111 @@ module sources
|
||||
dh(2) = ady(pdata%meta%level)
|
||||
dh(3) = adz(pdata%meta%level)
|
||||
|
||||
! add the EGLM-MHD source terms
|
||||
!
|
||||
if (glm_type > 0) then
|
||||
|
||||
! calculate the magnetic field divergence
|
||||
!
|
||||
call divergence(dh(:), pdata%q(ibx:ibz,:,:,:), db(:,:,:))
|
||||
|
||||
! update the momentum component increments, i.e.
|
||||
! d/dt (ρv) + ∇.F = - (∇.B)B
|
||||
!
|
||||
du(imx,:,:,:) = du(imx,:,:,:) - db(:,:,:) * pdata%q(ibx,:,:,:)
|
||||
du(imy,:,:,:) = du(imy,:,:,:) - db(:,:,:) * pdata%q(iby,:,:,:)
|
||||
du(imz,:,:,:) = du(imz,:,:,:) - db(:,:,:) * pdata%q(ibz,:,:,:)
|
||||
|
||||
! update the energy equation
|
||||
!
|
||||
if (ien > 0 .and. ibp > 0) then
|
||||
|
||||
! calculate the gradient of divergence potential
|
||||
!
|
||||
call gradient(dh(:), pdata%q(ibp,:,:,:), jc(inx:inz,:,:,:))
|
||||
|
||||
! add the divergence potential source term to the energy equation, i.e.
|
||||
! d/dt E + ∇.F = - B.(∇ψ)
|
||||
!
|
||||
du(ien,:,:,:) = du(ien,:,:,:) &
|
||||
- sum(pdata%q(ibx:ibz,:,:,:) * jc(inx:inz,:,:,:), 1)
|
||||
|
||||
end if ! ien > 0
|
||||
|
||||
! add the HEGLM-MHD source terms
|
||||
!
|
||||
if (glm_type > 1) then
|
||||
|
||||
! update magnetic field component increments, i.e.
|
||||
! d/dt B + ∇.F = - (∇.B)v
|
||||
!
|
||||
du(ibx,:,:,:) = du(ibx,:,:,:) - db(:,:,:) * pdata%q(ivx,:,:,:)
|
||||
du(iby,:,:,:) = du(iby,:,:,:) - db(:,:,:) * pdata%q(ivy,:,:,:)
|
||||
du(ibz,:,:,:) = du(ibz,:,:,:) - db(:,:,:) * pdata%q(ivz,:,:,:)
|
||||
|
||||
! update the energy equation
|
||||
!
|
||||
if (ien > 0) then
|
||||
|
||||
! calculate scalar product of velocity and magnetic field
|
||||
!
|
||||
jc(inx,:,:,:) = sum(pdata%q(ivx:ivz,:,:,:) &
|
||||
* pdata%q(ibx:ibz,:,:,:), 1)
|
||||
|
||||
! add the divergence potential source term to the energy equation, i.e.
|
||||
! d/dt E + ∇.F = - (∇.B) (v.B)
|
||||
!
|
||||
du(ien,:,:,:) = du(ien,:,:,:) - db(:,:,:) * jc(inx,:,:,:)
|
||||
|
||||
end if ! ien > 0
|
||||
|
||||
end if ! glm_type > 1
|
||||
|
||||
end if ! glmtype > 0
|
||||
|
||||
! proceed only if the resistivity coefficient is not zero
|
||||
!
|
||||
if (resistivity > 0.0d+00) then
|
||||
|
||||
! calculate the Laplace operator of B, i.e. Δ(B)
|
||||
!
|
||||
call laplace(dh(:), pdata%q(ibx,:,:,:), jc(inx,:,:,:))
|
||||
call laplace(dh(:), pdata%q(iby,:,:,:), jc(iny,:,:,:))
|
||||
call laplace(dh(:), pdata%q(ibz,:,:,:), jc(inz,:,:,:))
|
||||
call laplace(dh(:), pdata%q(ibx,:,:,:), jc(inx,:,:,:))
|
||||
call laplace(dh(:), pdata%q(iby,:,:,:), jc(iny,:,:,:))
|
||||
call laplace(dh(:), pdata%q(ibz,:,:,:), jc(inz,:,:,:))
|
||||
|
||||
! update magnetic field component increments
|
||||
!
|
||||
du(ibx,:,:,:) = du(ibx,:,:,:) + resistivity * jc(inx,:,:,:)
|
||||
du(iby,:,:,:) = du(iby,:,:,:) + resistivity * jc(iny,:,:,:)
|
||||
du(ibz,:,:,:) = du(ibz,:,:,:) + resistivity * jc(inz,:,:,:)
|
||||
du(ibx,:,:,:) = du(ibx,:,:,:) + resistivity * jc(inx,:,:,:)
|
||||
du(iby,:,:,:) = du(iby,:,:,:) + resistivity * jc(iny,:,:,:)
|
||||
du(ibz,:,:,:) = du(ibz,:,:,:) + resistivity * jc(inz,:,:,:)
|
||||
|
||||
! update energy equation
|
||||
!
|
||||
if (ien > 0) then
|
||||
if (ien > 0) then
|
||||
|
||||
! add the first resistive source term to the energy equation, i.e.
|
||||
! d/dt E + ∇.F = η B.[Δ(B)]
|
||||
!
|
||||
du(ien,:,:,:) = du(ien,:,:,:) &
|
||||
du(ien,:,:,:) = du(ien,:,:,:) &
|
||||
+ resistivity * (pdata%q(ibx,:,:,:) * jc(inx,:,:,:) &
|
||||
+ pdata%q(iby,:,:,:) * jc(iny,:,:,:) &
|
||||
+ pdata%q(ibz,:,:,:) * jc(inz,:,:,:))
|
||||
|
||||
! calculate current density J = ∇xB
|
||||
!
|
||||
call curl(dh(:), pdata%q(ibx:ibz,:,:,:), jc(inx:inz,:,:,:))
|
||||
call curl(dh(:), pdata%q(ibx:ibz,:,:,:), jc(inx:inz,:,:,:))
|
||||
|
||||
! add the second resistive source term to the energy equation, i.e.
|
||||
! d/dt E + ∇.F = η J²
|
||||
!
|
||||
du(ien,:,:,:) = du(ien,:,:,:) &
|
||||
du(ien,:,:,:) = du(ien,:,:,:) &
|
||||
+ resistivity * sum(jc(:,:,:,:) * jc(:,:,:,:), 1)
|
||||
|
||||
end if ! energy equation present
|
||||
end if ! energy equation present
|
||||
|
||||
end if ! resistivity is not zero
|
||||
end if ! resistivity is not zero
|
||||
|
||||
end if ! ibx > 0
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for source terms
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user