First step of implementation of time advance using new method.

BLOCK STRUCTURE

 - add new array to the data block structure to store the electromotive
   force components; these components are located at the centers of cell
   edges, in this way the CT update of the staggered magnetic field
   component will be easier;

EVOLUTION

 - add new subroutine advance() which performs several steps in order to
   advance the solution in time by one-step update; the substeps are the
   updates of the numerical flux, the flux boundary, the advance in time
   the solution, the updates of mesh structure, and the boundaries of
   conserved variables, and finally the new time step estimation;

 - add new subroutine update_flux() to update the numerical fluxes
   stored in the data blocks;

SCHEME

 - add new subroutine numerical_flux() to calculate fluxes at the proper
   locations;

 - add new logical argument to HLL and HLLC in order to specify if the
   flux should be returned as a numerical flux or its derivative;
This commit is contained in:
Grzegorz Kowal 2010-07-27 19:26:15 -03:00
parent 0b62ce8c35
commit ed88bc2526
3 changed files with 435 additions and 13 deletions

View File

@ -130,6 +130,7 @@ module blocks
real, dimension(:,:,:,:) , allocatable :: u ! the array of the conserved variables real, dimension(:,:,:,:) , allocatable :: u ! the array of the conserved variables
real, dimension(:,:,:,:,:), allocatable :: f ! the array of the numerical fluxes real, dimension(:,:,:,:,:), allocatable :: f ! the array of the numerical fluxes
real, dimension(:,:,:,:) , allocatable :: e ! the array of the electomagnetic force
end type block_data end type block_data
! define block_info structure for boundary exchange ! define block_info structure for boundary exchange
@ -370,6 +371,10 @@ module blocks
! !
allocate(pdata%f(NDIMS,nfl,im,jm,km)) allocate(pdata%f(NDIMS,nfl,im,jm,km))
! allocate space for the electromotive force
!
allocate(pdata%e( 3,im,jm,km))
! increase the number of allocated meta blocks ! increase the number of allocated meta blocks
! !
dblocks = dblocks + 1 dblocks = dblocks + 1
@ -504,6 +509,7 @@ module blocks
! deallocate fluxes ! deallocate fluxes
! !
deallocate(pdata%f) deallocate(pdata%f)
deallocate(pdata%e)
! nullify pointers ! nullify pointers
! !

View File

@ -35,6 +35,119 @@ module evolution
! !
!=============================================================================== !===============================================================================
! !
! advance: subroutine sweeps over all data blocks and updates their numerical
! fluxes; then updates the flux boundaries; next, advances in time the
! conserved variables, refines the mesh, updates the
! boundaries of conserved variables, and finally calculates the new
! time step
!
!===============================================================================
!
subroutine advance
use blocks , only : block_data, list_data
use boundaries , only : boundary!, boundary_flux
use mesh , only : update_mesh, dx_min
#ifdef MPI
use mpitools , only : mallreduceminr
#endif /* MPI */
use scheme , only : maxspeed!, advance_variables
use timer , only : start_timer, stop_timer
implicit none
! local variables
!
type(block_data), pointer :: pblock
real :: cmax, cm
!
!-------------------------------------------------------------------------------
!
! 1. iterate over all data blocks and calculate the numerical fluxes
!
pblock => list_data
do while (associated(pblock))
! - update the numerical fluxes of the current block
!
call update_flux(pblock)
! - assign pointer to the next block
!
pblock => pblock%next
end do
! ! 2. update the flux boundaries
! !
! call start_timer(4)
! call boundary_flux
! call stop_timer(4)
!
! ! 3. iterate over all data blocks and advance in time the conserved variables
! !
! pblock => list_data
! do while (associated(pblock))
!
! ! - advance in time the conserved variables of the current block
! !
! call advance_variables(pblock)
!
! ! - assign pointer to the next block
! !
! pblock => pblock%next
!
! end do
!
! ! 4. update the mesh (refine and derefine blocks)
! !
! call start_timer(5)
! call update_mesh
! call stop_timer(5)
!
! ! 5. update the boundaries of the conserved variables
! !
! call start_timer(4)
! call boundary
! call stop_timer(4)
!
! ! 6. iterate over all blocks in order to find the maximum speed
! !
! cmax = 1.0e-8
!
! pblock => list_data
! do while (associated(pblock))
!
! ! - calculate the maximum speed for the current block
! !
! cm = maxspeed(pblock%u)
!
! ! - take the maximum of the global and local maximum speeds
! !
! cmax = max(cmax, cm)
!
! ! - assign pointer to the next block
! !
! pblock => pblock%next
!
! end do
!
! ! - calculate the new time step
! !
! dtn = dx_min / max(cmax, 1.e-8)
!
! #ifdef MPI
! ! - reduce the new time step over all processes
! !
! call mallreduceminr(dtn)
! #endif /* MPI */
!
!-------------------------------------------------------------------------------
!
end subroutine advance
!
!===============================================================================
!
! evolve: subroutine sweeps over all leaf blocks and performs one step time ! evolve: subroutine sweeps over all leaf blocks and performs one step time
! evolution for each according to the selected integration scheme ! evolution for each according to the selected integration scheme
! !
@ -137,6 +250,35 @@ module evolution
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
end subroutine evolve end subroutine evolve
!
!===============================================================================
!
! update_flux: subroutine updates the fluxes for the current block using the
! conserved variables
!
!===============================================================================
!
subroutine update_flux(pblock)
use blocks , only : block_data, nqt
use config , only : im, jm, km
use scheme , only : numerical_flux
implicit none
! input arguments
!
type(block_data), intent(inout) :: pblock
!
!-------------------------------------------------------------------------------
!
! calculate the numerical flux for the block
!
call numerical_flux(pblock%u, pblock%f, pblock%e)
!
!-------------------------------------------------------------------------------
!
end subroutine update_flux
#ifdef EULER #ifdef EULER
! !
!=============================================================================== !===============================================================================

View File

@ -32,6 +32,270 @@ module scheme
! !
!============================================================================== !==============================================================================
! !
! numerical_flux: subroutine updates the numerical flux of the current data
! block
!
! notes: - the fluid fluxes must be located at cell interface centers
! [ fx(i+1/2,j,k), fy(i,j+1/2,k), fz(i,j,k+1/2) ]
! - the magnetic fluxes must be located at the cell edge centers
! [ ex(i,j+1/2,j+1/2), ey(i+1/2,j,j+1/2), ez(i+1/2,j+1/2,k) ]
!
!==============================================================================
!
subroutine numerical_flux(u, f, e)
use blocks , only : block_data
use blocks , only : nvr, nqt
use blocks , only : idn, imx, imy, imz
#ifdef ADI
use blocks , only : ien
#endif /* ADI */
#ifdef MHD
use blocks , only : ibx, iby, ibz, icx, icy, icz
#endif /* MHD */
use config , only : im, jm, km
#ifdef FLUXCT
use interpolation, only : magtocen
#endif /* FLUXCT */
implicit none
! input/output arguments
!
real, dimension( nqt,im,jm,km), intent(in) :: u
real, dimension(NDIMS,nqt,im,jm,km), intent(out) :: f
real, dimension( 3,im,jm,km), intent(out) :: e
! local variables
!
integer :: i, j, k
#ifdef MHD
integer :: im1, jm1, km1
#endif /* MHD */
! local temporary arrays
!
real, dimension(nvr,im) :: ux
real, dimension(nqt,im) :: fx
real, dimension(nvr,jm) :: uy
real, dimension(nqt,jm) :: fy
#if NDIMS == 3
real, dimension(nvr,km) :: uz
real, dimension(nqt,km) :: fz
#endif /* NDIMS == 3 */
real, dimension(nvr,im,jm,km) :: w
!
!-------------------------------------------------------------------------------
!
! 0. prepare the numerical flux calculation
!
! - reset the flux array to zero
!
f(:,:,:,:,:) = 0.0
e(:,:,:,:) = 0.0
! - copy the conserved variables to the local array
!
w(1:nqt,1:im,1:jm,1:km) = u(1:nqt,1:im,1:jm,1:km)
#ifdef FLUXCT
! - interpolate the cell centered components of the magnetic field
!
call magtocen(w)
#endif /* FLUXCT */
! 1. update along the X direction
!
do k = 1, km
#if NDIMS == 3 && defined MHD
km1 = max(k - 1, 1)
#endif /* NDIMS == 3 & MHD */
do j = 1, jm
#ifdef MHD
jm1 = max(j - 1, 1)
#endif /* MHD */
! - copy the directional vectors of variables for the one dimensional solver
!
ux(idn,1:im) = w(idn,1:im,j,k)
ux(imx,1:im) = w(imx,1:im,j,k)
ux(imy,1:im) = w(imy,1:im,j,k)
ux(imz,1:im) = w(imz,1:im,j,k)
#ifdef ADI
ux(ien,1:im) = w(ien,1:im,j,k)
#endif /* ADI */
#ifdef MHD
ux(ibx,1:im) = w(ibx,1:im,j,k)
ux(iby,1:im) = w(iby,1:im,j,k)
ux(ibz,1:im) = w(ibz,1:im,j,k)
#ifdef FLUXCT
ux(icx,1:im) = w(icx,1:im,j,k)
ux(icy,1:im) = w(icy,1:im,j,k)
ux(icz,1:im) = w(icz,1:im,j,k)
#endif /* FLUXCT */
#endif /* MHD */
! - execute the Riemann solver (returns fluxes for the update)
!
#ifdef HLL
call hll (im, ux, fx, .true.)
#endif /* HLL */
#ifdef HLLC
call hllc(im, ux, fx, .true.)
#endif /* HLLC */
! - update the fluxes along the current direction
!
f(1,idn,1:im,j,k) = fx(idn,1:im)
f(1,imx,1:im,j,k) = fx(imx,1:im)
f(1,imy,1:im,j,k) = fx(imy,1:im)
f(1,imz,1:im,j,k) = fx(imz,1:im)
#ifdef ADI
f(1,ien,1:im,j,k) = fx(ien,1:im)
#endif /* ADI */
#ifdef MHD
e(ibx,1:im,j ,k ) = e(ibx,1:im,j ,k ) + fx(ibx,1:im)
#if NDIMS == 3
e(iby,1:im,j ,k ) = e(iby,1:im,j ,k ) - 0.25 * fx(ibz,1:im)
e(iby,1:im,j ,km1) = e(iby,1:im,j ,km1) - 0.25 * fx(ibz,1:im)
#else /* NDIMS == 3 */
e(iby,1:im,j ,k ) = e(iby,1:im,j ,k ) - fx(ibz,1:im)
#endif /* NDIMS == 3 */
e(ibz,1:im,j ,k ) = e(ibz,1:im,j ,k ) + 0.25 * fx(iby,1:im)
e(ibz,1:im,jm1,k ) = e(ibz,1:im,jm1,k ) + 0.25 * fx(iby,1:im)
#endif /* MHD */
end do
end do
! 2. update along the Y direction
!
do k = 1, km
#if NDIMS == 3 && defined MHD
km1 = max(k - 1, 1)
#endif /* NDIMS == 3 & MHD */
do i = 1, im
#ifdef MHD
im1 = max(i - 1, 1)
#endif /* MHD */
! - copy the directional vectors of variables for the one dimensional solver
!
uy(idn,1:jm) = w(idn,i,1:jm,k)
uy(imx,1:jm) = w(imy,i,1:jm,k)
uy(imy,1:jm) = w(imz,i,1:jm,k)
uy(imz,1:jm) = w(imx,i,1:jm,k)
#ifdef ADI
uy(ien,1:jm) = w(ien,i,1:jm,k)
#endif /* ADI */
#ifdef MHD
uy(ibx,1:jm) = w(iby,i,1:jm,k)
uy(iby,1:jm) = w(ibz,i,1:jm,k)
uy(ibz,1:jm) = w(ibx,i,1:jm,k)
#ifdef FLUXCT
uy(icx,1:jm) = w(icy,i,1:jm,k)
uy(icy,1:jm) = w(icz,i,1:jm,k)
uy(icz,1:jm) = w(icx,i,1:jm,k)
#endif /* FLUXCT */
#endif /* MHD */
! - execute the Riemann solver (returns fluxes for the update)
!
#ifdef HLL
call hll (jm, uy, fy, .true.)
#endif /* HLL */
#ifdef HLLC
call hllc(jm, uy, fy, .true.)
#endif /* HLLC */
! - update the fluxes along the current direction
!
f(2,idn,i,1:jm,k) = fy(idn,1:jm)
f(2,imx,i,1:jm,k) = fy(imz,1:jm)
f(2,imy,i,1:jm,k) = fy(imx,1:jm)
f(2,imz,i,1:jm,k) = fy(imy,1:jm)
#ifdef ADI
f(2,ien,i,1:jm,k) = fy(ien,1:jm)
#endif /* ADI */
#ifdef MHD
#if NDIMS == 3
e(ibx,i ,1:jm,k ) = e(ibx,i ,1:jm,k ) + 0.25 * fy(ibz,1:jm)
e(ibx,i ,1:jm,km1) = e(ibx,i ,1:jm,km1) + 0.25 * fy(ibz,1:jm)
#else /* NDIMS == 3 */
e(ibx,i ,1:jm,k ) = e(ibx,i ,1:jm,k ) + fy(ibz,1:jm)
#endif /* NDIMS == 3 */
e(iby,i ,1:jm,k ) = e(iby,i ,1:jm,k ) + fy(iby,1:jm)
e(ibz,i ,1:jm,k ) = e(ibz,i ,1:jm,k ) - 0.25 * fy(ibx,1:jm)
e(ibz,im1,1:jm,k ) = e(ibz,im1,1:jm,k ) - 0.25 * fy(ibx,1:jm)
#endif /* MHD */
end do
end do
#if NDIMS == 3
! 3. update along the Z direction
!
do j = 1, jm
#ifdef MHD
jm1 = max(j - 1, 1)
#endif /* MHD */
do i = 1, im
#ifdef MHD
im1 = max(i - 1, 1)
#endif /* MHD */
! - copy the directional vectors of variables for the one dimensional solver
!
uz(idn,1:km) = w(idn,i,j,1:km)
uz(imx,1:km) = w(imz,i,j,1:km)
uz(imy,1:km) = w(imx,i,j,1:km)
uz(imz,1:km) = w(imy,i,j,1:km)
#ifdef ADI
uz(ien,1:km) = w(ien,i,j,1:km)
#endif /* ADI */
#ifdef MHD
uz(ibx,1:km) = w(ibz,i,j,1:km)
uz(iby,1:km) = w(ibx,i,j,1:km)
uz(ibz,1:km) = w(iby,i,j,1:km)
#ifdef FLUXCT
uz(icx,1:km) = w(icz,i,j,1:km)
uz(icy,1:km) = w(icx,i,j,1:km)
uz(icz,1:km) = w(icy,i,j,1:km)
#endif /* FLUXCT */
#endif /* MHD */
! - execute the Riemann solver (returns fluxes for the update)
!
#ifdef HLL
call hll (km, uz, fz, .true.)
#endif /* HLL */
#ifdef HLLC
call hllc(km, uz, fz, .true.)
#endif /* HLLC */
! - update the fluxes along the current direction
!
f(3,idn,i,j,1:km) = fz(idn,1:km)
f(3,imx,i,j,1:km) = fz(imy,1:km)
f(3,imy,i,j,1:km) = fz(imz,1:km)
f(3,imz,i,j,1:km) = fz(imx,1:km)
#ifdef ADI
f(3,ien,i,j,1:km) = fz(ien,1:km)
#endif /* ADI */
#ifdef MHD
e(ibx,i ,j ,1:km) = e(ibx,i ,j ,1:km) - 0.25 * fz(iby,1:km)
e(ibx,i ,jm1,1:km) = e(ibx,i ,jm1,1:km) - 0.25 * fz(iby,1:km)
e(iby,i ,j ,1:km) = e(iby,i ,j ,1:km) + 0.25 * fz(ibx,1:km)
e(iby,im1,j ,1:km) = e(iby,im1,j ,1:km) + 0.25 * fz(ibx,1:km)
e(ibz,i ,j ,1:km) = e(ibz,i ,j ,1:km) + fz(ibz,1:km)
#endif /* MHD */
end do
end do
#endif /* NDIMS == 3 */
!
!-------------------------------------------------------------------------------
!
end subroutine numerical_flux
!
!==============================================================================
!
! update: subroutine sweeps over all directions and integrates the directional ! update: subroutine sweeps over all directions and integrates the directional
! derivatives of the flux in order to get the increment of solution ! derivatives of the flux in order to get the increment of solution
! !
@ -57,7 +321,7 @@ module scheme
! input arguments ! input arguments
! !
real, dimension(nqt,im,jm,km) , intent(in) :: u real, dimension(nqt,im,jm,km) , intent(in) :: u
real, dimension(NDIMS,nfl,im,jm,km), intent(out) :: f real, dimension(NDIMS,nqt,im,jm,km), intent(out) :: f
real, dimension(nqt,im,jm,km) , intent(out) :: du real, dimension(nqt,im,jm,km) , intent(out) :: du
real , intent(in) :: dxi, dyi, dzi real , intent(in) :: dxi, dyi, dzi
@ -132,10 +396,10 @@ module scheme
! execute solver (returns fluxes for the update) ! execute solver (returns fluxes for the update)
! !
#ifdef HLL #ifdef HLL
call hll (im, ux, fx) call hll (im, ux, fx, .false.)
#endif /* HLL */ #endif /* HLL */
#ifdef HLLC #ifdef HLLC
call hllc(im, ux, fx) call hllc(im, ux, fx, .false.)
#endif /* HLLC */ #endif /* HLLC */
! update fluxes along the X direction ! update fluxes along the X direction
@ -232,10 +496,10 @@ module scheme
! execute solver (returns fluxes for the update) ! execute solver (returns fluxes for the update)
! !
#ifdef HLL #ifdef HLL
call hll (jm, uy, fy) call hll (jm, uy, fy, .false.)
#endif /* HLL */ #endif /* HLL */
#ifdef HLLC #ifdef HLLC
call hllc(jm, uy, fy) call hllc(jm, uy, fy, .false.)
#endif /* HLLC */ #endif /* HLLC */
! update fluxes along the Y direction ! update fluxes along the Y direction
@ -331,10 +595,10 @@ module scheme
! execute solver (returns fluxes for the update) ! execute solver (returns fluxes for the update)
! !
#ifdef HLL #ifdef HLL
call hll (km, uz, fz) call hll (km, uz, fz, .false.)
#endif /* HLL */ #endif /* HLL */
#ifdef HLLC #ifdef HLLC
call hllc(km, uz, fz) call hllc(km, uz, fz, .false.)
#endif /* HLLC */ #endif /* HLLC */
! update fluxes along the Z direction ! update fluxes along the Z direction
@ -400,7 +664,7 @@ module scheme
! !
!=============================================================================== !===============================================================================
! !
subroutine hll(n, u, f) subroutine hll(n, u, f, d)
use blocks , only : nvr, nfl, nqt use blocks , only : nvr, nfl, nqt
use blocks , only : ivx, ivz use blocks , only : ivx, ivz
@ -416,6 +680,7 @@ module scheme
integer , intent(in) :: n integer , intent(in) :: n
real, dimension(nvr,n), intent(in) :: u real, dimension(nvr,n), intent(in) :: u
real, dimension(nqt,n), intent(out) :: f real, dimension(nqt,n), intent(out) :: f
logical , intent(in) :: d
! local variables ! local variables
! !
@ -486,15 +751,19 @@ module scheme
! calculate numerical flux ! calculate numerical flux
! !
f( 1:nfl,2:n) = - fn( 1:nfl,2:n) + fn( 1:nfl,1:n-1) if (d) then
f(1:nqt,1:n) = fn(1:nqt,1:n)
else
f( 1:nfl,2:n) = - fn( 1:nfl,2:n) + fn( 1:nfl,1:n-1)
#ifdef MHD #ifdef MHD
#ifdef FIELDCD #ifdef FIELDCD
call emf(n, q(ivx:ivz,:), q(ibx:ibz,:), f(ibx:ibz,:)) call emf(n, q(ivx:ivz,:), q(ibx:ibz,:), f(ibx:ibz,:))
#endif /* FIELDCD */ #endif /* FIELDCD */
#ifdef FLUXCT #ifdef FLUXCT
f(ibx:ibz,1:n) = 0.25 * fn(ibx:ibz,1:n) f(ibx:ibz,1:n) = 0.25 * fn(ibx:ibz,1:n)
#endif /* FLUXCT */ #endif /* FLUXCT */
#endif /* MHD */ #endif /* MHD */
end if
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
@ -509,7 +778,7 @@ module scheme
! !
!=============================================================================== !===============================================================================
! !
subroutine hllc(n, u, f) subroutine hllc(n, u, f, d)
use blocks , only : nvr, nfl, nqt use blocks , only : nvr, nfl, nqt
use interpolation, only : reconstruct use interpolation, only : reconstruct
@ -521,6 +790,7 @@ module scheme
integer , intent(in) :: n integer , intent(in) :: n
real, dimension(nvr,n), intent(in) :: u real, dimension(nvr,n), intent(in) :: u
real, dimension(nqt,n), intent(out) :: f real, dimension(nqt,n), intent(out) :: f
logical , intent(in) :: d
! local variables ! local variables
! !
@ -691,7 +961,11 @@ module scheme
! calculate numerical flux ! calculate numerical flux
! !
f(:,2:n) = - fx(:,2:n) + fx(:,1:n-1) if (d) then
f(:,:) = fx(:,:)
else
f(:,2:n) = - fx(:,2:n) + fx(:,1:n-1)
end if
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !