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

1180 lines
31 KiB
Fortran
Raw Normal View History

Update the copyright info and fix the length of separating lines.
2010-10-13 03:32:10 -03:00
!!******************************************************************************
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
!!
Update copyright headers with the new email.
2011-05-05 18:37:53 -03:00
!! module: EVOLUTION - handling the time evolution of the block structure
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
!!
Update copyright headers with the new email.
2011-05-05 18:37:53 -03:00
!! Copyright (C) 2008-2011 Grzegorz Kowal <grzegorz@amuncode.org>
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
!!
Update the copyright info and fix the length of separating lines.
2010-10-13 03:32:10 -03:00
!!******************************************************************************
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
!!
Update copyright headers with the new email.
2011-05-05 18:37:53 -03:00
!! This file is part of the AMUN code.
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
!!
Change the license to GNU General Public License, version 2.
2011-04-29 11:21:30 -03:00
!! This program is free software; you can redistribute it and/or
!! modify it under the terms of the GNU General Public License
!! as published by the Free Software Foundation; either version 2
!! of the License, or (at your option) any later version.
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
!!
Change the license to GNU General Public License, version 2.
2011-04-29 11:21:30 -03:00
!! This program is distributed in the hope that it will be useful,
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
Change the license to GNU General Public License, version 2.
2011-04-29 11:21:30 -03:00
!! along with this program; if not, write to the Free Software Foundation,
!! Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
!!
Update the copyright info and fix the length of separating lines.
2010-10-13 03:32:10 -03:00
!!******************************************************************************
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
!!
!
module evolution
implicit none
integer, save :: n
SCHEME: move cmax parameter from module 'evolution' to 'scheme'.
2010-12-01 13:13:27 -02:00
real , save :: t, dt, dtn
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
contains
!
Updated comments in the evolution module.
2008-12-08 21:07:10 -06:00
!===============================================================================
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
!
Updated comments in the evolution module.
2008-12-08 21:07:10 -06:00
! evolve: subroutine sweeps over all leaf blocks and performs one step time
! evolution for each according to the selected integration scheme
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
!
Updated comments in the evolution module.
2008-12-08 21:07:10 -06:00
!===============================================================================
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
!
EVOLUTION: remove unused subroutines. - remove subroutine advance(); this subroutine is good for GALERKIN approach; it's gonna be reintroduced during this method implementation; - remove subroutine update_flux(); this subroutine is called from advance() and is not used in the GLM-MHD approach; - make GNU Fortran happy bu adding () to the subroutine calls;
2010-12-01 15:14:07 -02:00
subroutine evolve()
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
EVOLUTION: evolve the forcing source terms. - call subroutine evolve_forcing() before the update of all blocks; this subroutine evolves the forcing source terms by an interval dt in the Fourier space; then during the update the forcing Fourier coefficients will be transformed to real space for each block separately;
2010-12-08 18:19:16 -02:00
use blocks , only : block_data, list_data
use boundaries, only : boundary_variables
Implement subroutines to correct boundary fluxes. - subroutine boundary_correct_fluxes() sweeps over all leafs and if it finds that a neighbor is at higher level, performs a correction of the numerical flux of current block in this direction; this subroutine does not support MPI yet!; - subroutine correct_flux() performs the atual correction of the numerical flux; only part belonging to domain is updated, since ghost zones are updated anyway in the next step during the variable update; - call boundary_correct_flux() from evolve();
2011-04-30 15:10:02 -03:00
#ifdef CONSERVATIVE
use boundaries, only : boundary_correct_fluxes
#endif /* CONSERVATIVE */
Perform the refinement step only if maxlev > 1.
2011-04-26 12:37:25 -03:00
#ifdef REFINE
use config , only : maxlev
#endif /* REFINE */
Implement viscous terms.
2011-03-18 15:33:24 -03:00
#ifdef VISCOSITY
use config , only : visc
#endif /* VISCOSITY */
Change flag RESIS to RESISTIVITY.
2011-03-22 17:24:28 -03:00
#if defined MHD && defined RESISTIVITY
EVOLUTION: take resistivity timestep condition into account.
2010-12-14 17:29:38 -02:00
use config , only : ueta
Change flag RESIS to RESISTIVITY.
2011-03-22 17:24:28 -03:00
#endif /* MHD & RESISTIVITY */
EVOLUTION: evolve the forcing source terms. - call subroutine evolve_forcing() before the update of all blocks; this subroutine evolves the forcing source terms by an interval dt in the Fourier space; then during the update the forcing Fourier coefficients will be transformed to real space for each block separately;
2010-12-08 18:19:16 -02:00
#ifdef FORCE
FORCE: implement and use Fourier transform of velocity. - we need to calculate Fourier transform of velocity for all forcing components in order to minimize the velocity-force correlation;
2011-03-07 00:08:31 -03:00
use forcing , only : fourier_transform, evolve_forcing
EVOLUTION: evolve the forcing source terms. - call subroutine evolve_forcing() before the update of all blocks; this subroutine evolves the forcing source terms by an interval dt in the Fourier space; then during the update the forcing Fourier coefficients will be transformed to real space for each block separately;
2010-12-08 18:19:16 -02:00
#endif /* FORCE */
use mesh , only : update_mesh
use mesh , only : dx_min
use scheme , only : cmax
Start/stop timers from the inside of evolution subroutines.
2011-05-03 00:01:00 -03:00
use timer , only : start_timer, stop_timer
FORCE: implement and use Fourier transform of velocity. - we need to calculate Fourier transform of velocity for all forcing components in order to minimize the velocity-force correlation;
2011-03-07 00:08:31 -03:00
use variables , only : idn, imz
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
implicit none
! local variables
!
Evolve new list of data blocks. TIME INTEGRATION - now update the solution using new list of data blocks belonging to the current process only
2009-09-14 18:28:17 -03:00
type(block_data), pointer :: pblock
EVOLUTION: add subroutine to update the maximum speed. - make the variable cmax global in the module 'evolution'; - add a new subroutine update_maximum_speed() which updates the maximum speed cmax in the system iterating over all data blocks; - use the subroutine update_maximum_speed() in evolve();
2010-12-01 10:39:18 -02:00
real :: cm
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
!
Updated comments in the evolution module.
2008-12-08 21:07:10 -06:00
!-------------------------------------------------------------------------------
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
!
Start/stop timers from the inside of evolution subroutines.
2011-05-03 00:01:00 -03:00
! start the evolution timer
!
call start_timer(2)
EVOLUTION: evolve the forcing source terms. - call subroutine evolve_forcing() before the update of all blocks; this subroutine evolves the forcing source terms by an interval dt in the Fourier space; then during the update the forcing Fourier coefficients will be transformed to real space for each block separately;
2010-12-08 18:19:16 -02:00
#ifdef FORCE
FORCE: implement and use Fourier transform of velocity. - we need to calculate Fourier transform of velocity for all forcing components in order to minimize the velocity-force correlation;
2011-03-07 00:08:31 -03:00
! perform the Fourier transform of the velocity field
!
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,:,:,:))
pblock => pblock%next ! assign pointer to the next block
end do
EVOLUTION: evolve the forcing source terms. - call subroutine evolve_forcing() before the update of all blocks; this subroutine evolves the forcing source terms by an interval dt in the Fourier space; then during the update the forcing Fourier coefficients will be transformed to real space for each block separately;
2010-12-08 18:19:16 -02:00
! evolve the forcing source terms by the time interval dt
!
call evolve_forcing(dt)
#endif /* FORCE */
Evolve new list of data blocks. TIME INTEGRATION - now update the solution using new list of data blocks belonging to the current process only
2009-09-14 18:28:17 -03:00
! iterate over all data blocks and perform one step of time evolution
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
!
Evolve new list of data blocks. TIME INTEGRATION - now update the solution using new list of data blocks belonging to the current process only
2009-09-14 18:28:17 -03:00
pblock => list_data
Added CFL condition and calculation of the new time step. The function to calculate the maximum speed in the block has been added. This function is used to determine the maximum speed globally, which is next used to estimate the next time step.
2008-12-09 14:51:33 -06:00
do while (associated(pblock))
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
! check if this block is a leaf
!
Add a subroutine to calculate numerical fluxes. - the subroutine update_flux() in the SCHEME module calculates numerical fluxes from the conserved variables; the calculated fluxes are then returned;
2011-04-30 12:28:02 -03:00
#ifdef CONSERVATIVE
Implement flux integration using Euler method. - subroutine flux_euler() integrates the numerical flux using the first order Euler method;
2011-04-30 12:53:17 -03:00
if (pblock%meta%leaf) &
#ifdef EULER
call flux_euler(pblock)
#endif /* EULER */
Implement flux integration using RK2. - implement subroutine flux_rk2() which calculates fluxes at time t and t + dt and averages them; then the calculated flux is passed to variable update; - implement subroutine advance_solution() which advances an array of conserved variables using fluxes provided as an argument; - use advance_solution() in update_solution();
2011-05-01 08:57:14 -03:00
#ifdef RK2
call flux_rk2(pblock)
#endif /* RK2 */
Implement flux integration using RK3 method.
2011-05-01 10:31:27 -03:00
#ifdef RK3
call flux_rk3(pblock)
#endif /* RK3 */
Add a subroutine to calculate numerical fluxes. - the subroutine update_flux() in the SCHEME module calculates numerical fluxes from the conserved variables; the calculated fluxes are then returned;
2011-04-30 12:28:02 -03:00
#else /* CONSERVATIVE */
Evolve new list of data blocks. TIME INTEGRATION - now update the solution using new list of data blocks belonging to the current process only
2009-09-14 18:28:17 -03:00
if (pblock%meta%leaf) &
New time integration method EULER.
2008-12-09 21:06:21 -06:00
#ifdef EULER
call evolve_euler(pblock)
#endif /* EULER */
Runge-Kutta 2nd order time integration implemented. In addition, I've done some fixes to the problem initialization, and I defined new variables igrids, jgrids, kgrids, which specify the dimensions of the block.
2008-12-08 16:21:59 -06:00
#ifdef RK2
Added CFL condition and calculation of the new time step. The function to calculate the maximum speed in the block has been added. This function is used to determine the maximum speed globally, which is next used to estimate the next time step.
2008-12-09 14:51:33 -06:00
call evolve_rk2(pblock)
Runge-Kutta 2nd order time integration implemented. In addition, I've done some fixes to the problem initialization, and I defined new variables igrids, jgrids, kgrids, which specify the dimensions of the block.
2008-12-08 16:21:59 -06:00
#endif /* RK2 */
EVOLUTION: implement the 3rd order RK method. - implement the 3rd order Runge-Kutta time integration; - rewrite slightly the 2nd order Runge-Kutta method;
2010-12-03 15:38:48 -02:00
#ifdef RK3
call evolve_rk3(pblock)
#endif /* RK3 */
Add a subroutine to calculate numerical fluxes. - the subroutine update_flux() in the SCHEME module calculates numerical fluxes from the conserved variables; the calculated fluxes are then returned;
2011-04-30 12:28:02 -03:00
#endif /* CONSERVATIVE */
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
! assign pointer to the next block
!
Added CFL condition and calculation of the new time step. The function to calculate the maximum speed in the block has been added. This function is used to determine the maximum speed globally, which is next used to estimate the next time step.
2008-12-09 14:51:33 -06:00
pblock => pblock%next
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
Added CFL condition and calculation of the new time step. The function to calculate the maximum speed in the block has been added. This function is used to determine the maximum speed globally, which is next used to estimate the next time step.
2008-12-09 14:51:33 -06:00
end do
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
Add one step update of conserved variables. - after calculating an integrals of fluxes with a given order we performed one step updated of the conserved variables stored in data blocks using corresponding numerical flux stored in the same data block; the subroutine which performs update is called update_solution();
2011-04-30 13:14:07 -03:00
#ifdef CONSERVATIVE
Implement subroutines to correct boundary fluxes. - subroutine boundary_correct_fluxes() sweeps over all leafs and if it finds that a neighbor is at higher level, performs a correction of the numerical flux of current block in this direction; this subroutine does not support MPI yet!; - subroutine correct_flux() performs the atual correction of the numerical flux; only part belonging to domain is updated, since ghost zones are updated anyway in the next step during the variable update; - call boundary_correct_flux() from evolve();
2011-04-30 15:10:02 -03:00
! correct the numerical fluxes between neighboring blocks which are at different
! levels
!
Put condition from previous commit in boundary_correct_fluxes().
2011-05-01 10:44:31 -03:00
call boundary_correct_fluxes()
Implement subroutines to correct boundary fluxes. - subroutine boundary_correct_fluxes() sweeps over all leafs and if it finds that a neighbor is at higher level, performs a correction of the numerical flux of current block in this direction; this subroutine does not support MPI yet!; - subroutine correct_flux() performs the atual correction of the numerical flux; only part belonging to domain is updated, since ghost zones are updated anyway in the next step during the variable update; - call boundary_correct_flux() from evolve();
2011-04-30 15:10:02 -03:00
Add one step update of conserved variables. - after calculating an integrals of fluxes with a given order we performed one step updated of the conserved variables stored in data blocks using corresponding numerical flux stored in the same data block; the subroutine which performs update is called update_solution();
2011-04-30 13:14:07 -03:00
! update solution using numerical fluxes stored in data blocks
!
pblock => list_data
do while (associated(pblock))
! check if this block is a leaf and update its conserved variables using
! corrected numerical fluxes
!
if (pblock%meta%leaf) call update_solution(pblock)
! assign pointer to the next block
!
pblock => pblock%next
end do
#endif /* CONSERVATIVE */
New module to handle the boundary conditions. The subroutine 'boundary' sweeps over all leaf blocks. For each block it sweeps over its neighbors and performs update of the boundaries. This is an initial version yet, it supports only neighboring blocks of the same level of refinement.
2008-12-09 20:37:31 -06:00
! update boundaries
Store all variables in one array U. One array U, which is a field in the BLOCK structure, stores all variables. The number of variables is determined by the nvars parameter. To access each variable we use variable index now, like idn, imx, imy, mz, ien, etc.
2008-12-08 12:14:13 -06:00
!
EVOLUTION: remove unused subroutines. - remove subroutine advance(); this subroutine is good for GALERKIN approach; it's gonna be reintroduced during this method implementation; - remove subroutine update_flux(); this subroutine is called from advance() and is not used in the GLM-MHD approach; - make GNU Fortran happy bu adding () to the subroutine calls;
2010-12-01 15:14:07 -02:00
call boundary_variables()
Store all variables in one array U. One array U, which is a field in the BLOCK structure, stores all variables. The number of variables is determined by the nvars parameter. To access each variable we use variable index now, like idn, imx, imy, mz, ien, etc.
2008-12-08 12:14:13 -06:00
Introduce a compilation flag REFINE. - the flag REFINE determines if the code should be compiler with the adaptive mesh refining/derefining; otherwise, the static mesh is generated;
2011-03-10 01:08:00 -03:00
#ifdef REFINE
Perform the refinement step only if maxlev > 1.
2011-04-26 12:37:25 -03:00
! chec if we need to perform the refinement step
!
if (maxlev .gt. 1) then
Use TVD interpolation in boundary update, update mesh before time step. BOUNDARY CONDITIONS - use TVD interpolation for prolongation of the boundary conditions CONFIGURATION - put lower limit for the number of ghost and domain cells HOST FILES - host files should not be included in the revision control INTERPOLATION - define all arrays as REAL, not REAL(KIND=8) since the precision of calculations is determined at the compilation stage - replace j0 and j1 indices with new more obvious il and ir MESH - update mesh before calculating new time step
2010-03-14 15:40:24 -03:00
! check refinement and refine
!
Perform the refinement step only if maxlev > 1.
2011-04-26 12:37:25 -03:00
call update_mesh()
Use TVD interpolation in boundary update, update mesh before time step. BOUNDARY CONDITIONS - use TVD interpolation for prolongation of the boundary conditions CONFIGURATION - put lower limit for the number of ghost and domain cells HOST FILES - host files should not be included in the revision control INTERPOLATION - define all arrays as REAL, not REAL(KIND=8) since the precision of calculations is determined at the compilation stage - replace j0 and j1 indices with new more obvious il and ir MESH - update mesh before calculating new time step
2010-03-14 15:40:24 -03:00
! update boundaries
!
Perform the refinement step only if maxlev > 1.
2011-04-26 12:37:25 -03:00
call boundary_variables()
end if ! maxlev > 1
Introduce a compilation flag REFINE. - the flag REFINE determines if the code should be compiler with the adaptive mesh refining/derefining; otherwise, the static mesh is generated;
2011-03-10 01:08:00 -03:00
#endif /* REFINE */
Start/stop timers from the inside of evolution subroutines.
2011-05-03 00:01:00 -03:00
EVOLUTION: add subroutine to update the maximum speed. - make the variable cmax global in the module 'evolution'; - add a new subroutine update_maximum_speed() which updates the maximum speed cmax in the system iterating over all data blocks; - use the subroutine update_maximum_speed() in evolve();
2010-12-01 10:39:18 -02:00
! update the maximum speed
Added CFL condition and calculation of the new time step. The function to calculate the maximum speed in the block has been added. This function is used to determine the maximum speed globally, which is next used to estimate the next time step.
2008-12-09 14:51:33 -06:00
!
EVOLUTION: add subroutine to update the maximum speed. - make the variable cmax global in the module 'evolution'; - add a new subroutine update_maximum_speed() which updates the maximum speed cmax in the system iterating over all data blocks; - use the subroutine update_maximum_speed() in evolve();
2010-12-01 10:39:18 -02:00
call update_maximum_speed()
! get maximum time step
!
dtn = dx_min / max(cmax, 1.0d-16)
Implement viscous terms.
2011-03-18 15:33:24 -03:00
#ifdef VISCOSITY
dtn = min(dtn, 0.5d0 * dx_min * dx_min / max(1.0d-16, visc))
#endif /* VISCOSITY */
Change flag RESIS to RESISTIVITY.
2011-03-22 17:24:28 -03:00
#if defined MHD && defined RESISTIVITY
EVOLUTION: take resistivity timestep condition into account.
2010-12-14 17:29:38 -02:00
dtn = min(dtn, 0.5d0 * dx_min * dx_min / max(1.0d-16, ueta))
Change flag RESIS to RESISTIVITY.
2011-03-22 17:24:28 -03:00
#endif /* MHD & RESISTIVITY */
Start/stop timers from the inside of evolution subroutines.
2011-05-03 00:01:00 -03:00
! stop the evolution timer
!
call stop_timer(2)
EVOLUTION: add subroutine to update the maximum speed. - make the variable cmax global in the module 'evolution'; - add a new subroutine update_maximum_speed() which updates the maximum speed cmax in the system iterating over all data blocks; - use the subroutine update_maximum_speed() in evolve();
2010-12-01 10:39:18 -02:00
!
!-------------------------------------------------------------------------------
!
end subroutine evolve
!
!===============================================================================
!
! update_maximum_speed: subroutine updates module variable cmax with the value
! corresponding to the maximum speed in the system
!
!===============================================================================
!
subroutine update_maximum_speed()
EVOLUTION: reduce cmax instead of dt. - now we reduce the maximum speed over all processors instead of the time step;
2010-12-03 18:02:07 -02:00
use blocks , only : block_data, list_data
#ifdef MPI
use mpitools, only : mallreducemaxr
#endif /* MPI */
use scheme , only : maxspeed, cmax
Start/stop timers from the inside of evolution subroutines.
2011-05-03 00:01:00 -03:00
use timer , only : start_timer, stop_timer
EVOLUTION: add subroutine to update the maximum speed. - make the variable cmax global in the module 'evolution'; - add a new subroutine update_maximum_speed() which updates the maximum speed cmax in the system iterating over all data blocks; - use the subroutine update_maximum_speed() in evolve();
2010-12-01 10:39:18 -02:00
implicit none
! local variables
!
type(block_data), pointer :: pblock
real :: cm
!
!-------------------------------------------------------------------------------
!
Start/stop timers from the inside of evolution subroutines.
2011-05-03 00:01:00 -03:00
! start the evolution timer
!
call start_timer(6)
EVOLUTION: add subroutine to update the maximum speed. - make the variable cmax global in the module 'evolution'; - add a new subroutine update_maximum_speed() which updates the maximum speed cmax in the system iterating over all data blocks; - use the subroutine update_maximum_speed() in evolve();
2010-12-01 10:39:18 -02:00
! reset the maximum speed
!
cmax = 1.0d-16
Added CFL condition and calculation of the new time step. The function to calculate the maximum speed in the block has been added. This function is used to determine the maximum speed globally, which is next used to estimate the next time step.
2008-12-09 14:51:33 -06:00
! iterate over all blocks in order to find the maximum speed
!
Evolve new list of data blocks. TIME INTEGRATION - now update the solution using new list of data blocks belonging to the current process only
2009-09-14 18:28:17 -03:00
pblock => list_data
Added CFL condition and calculation of the new time step. The function to calculate the maximum speed in the block has been added. This function is used to determine the maximum speed globally, which is next used to estimate the next time step.
2008-12-09 14:51:33 -06:00
do while (associated(pblock))
! check if this block is a leaf
!
Evolve new list of data blocks. TIME INTEGRATION - now update the solution using new list of data blocks belonging to the current process only
2009-09-14 18:28:17 -03:00
if (pblock%meta%leaf) &
Added CFL condition and calculation of the new time step. The function to calculate the maximum speed in the block has been added. This function is used to determine the maximum speed globally, which is next used to estimate the next time step.
2008-12-09 14:51:33 -06:00
cm = maxspeed(pblock%u)
! compare global and local maximum speeds
!
cmax = max(cmax, cm)
! assign pointer to the next block
!
pblock => pblock%next
end do
EVOLUTION: reduce cmax instead of dt. - now we reduce the maximum speed over all processors instead of the time step;
2010-12-03 18:02:07 -02:00
#ifdef MPI
! reduce the maximum speed over all processes
!
call mallreducemaxr(cmax)
#endif /* MPI */
Start/stop timers from the inside of evolution subroutines.
2011-05-03 00:01:00 -03:00
! stop the evolution timer
!
call stop_timer(6)
Implement FLUXCT integration of the induction equation. SCHEME - implement Flux-CT scheme for the staggered magnetic field integration; BLOCK STRUCTURE - use more space efficient storage of the variables, which means storing only staggered components of magnetic field; cell-centered components are calculated only when necessary; EVOLUTION - remove loops in the field updates; operations are performed on the arrays; BOUNDARY CONDITIONS - remove loops in the bnd_copy; operations are calculated on the whole array now; INTERPOLATION - subroutine magtocen() has been rewritten to avoid problems with the array allocation; now as an argument we enter the array of all variables; subroutine uses indices for the face-centered and cell-centered magnetic field components internally; MAKE - add flag defining Flux-CT scheme; PROBLEM - use predefined array variables instead of allocated;
2010-02-28 18:35:57 -03:00
!
Updated comments in the evolution module.
2008-12-08 21:07:10 -06:00
!-------------------------------------------------------------------------------
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
!
EVOLUTION: add subroutine to update the maximum speed. - make the variable cmax global in the module 'evolution'; - add a new subroutine update_maximum_speed() which updates the maximum speed cmax in the system iterating over all data blocks; - use the subroutine update_maximum_speed() in evolve();
2010-12-01 10:39:18 -02:00
end subroutine update_maximum_speed
Add a subroutine to calculate numerical fluxes. - the subroutine update_flux() in the SCHEME module calculates numerical fluxes from the conserved variables; the calculated fluxes are then returned;
2011-04-30 12:28:02 -03:00
#ifdef CONSERVATIVE
Add one step update of conserved variables. - after calculating an integrals of fluxes with a given order we performed one step updated of the conserved variables stored in data blocks using corresponding numerical flux stored in the same data block; the subroutine which performs update is called update_solution();
2011-04-30 13:14:07 -03:00
!
!===============================================================================
!
! update_solution: subroutine performs an one step update of the conserved
! variables for the given data block using the integrated
! numerical fluxes stored in the same data block
!
!===============================================================================
!
subroutine update_solution(pblock)
use blocks , only : block_data
use config , only : im, jm, km
use mesh , only : adxi, adyi, adzi
Add magnetic terms to advance_solution() and update_solution().
2011-05-01 09:44:19 -03:00
#if defined MHD && defined GLM
use config , only : alpha_p
use mesh , only : dx_min
use scheme , only : cmax
use variables, only : iph
#endif /* MHD & GLM */
Include forcing term in update_solution().
2011-05-01 09:24:11 -03:00
#ifdef FORCE
Add magnetic terms to advance_solution() and update_solution().
2011-05-01 09:44:19 -03:00
use forcing , only : real_forcing
Include forcing term in update_solution().
2011-05-01 09:24:11 -03:00
use variables, only : inx, iny, inz
use variables, only : idn, imx, imy, imz
Add injected energy to total one in adiabatic case. - if we use forcing in the case of adiabatic equation of state, the total energy of each cell must be updated by the kinetic energy injected to this cell; otherwise the energy conservation will be violated;
2011-05-01 19:21:43 -03:00
#ifdef ADI
use variables, only : ien
#endif /* ADI */
Include forcing term in update_solution().
2011-05-01 09:24:11 -03:00
#endif /* FORCE */
Add one step update of conserved variables. - after calculating an integrals of fluxes with a given order we performed one step updated of the conserved variables stored in data blocks using corresponding numerical flux stored in the same data block; the subroutine which performs update is called update_solution();
2011-04-30 13:14:07 -03:00
implicit none
! input arguments
!
type(block_data), intent(inout) :: pblock
Implement flux integration using RK2. - implement subroutine flux_rk2() which calculates fluxes at time t and t + dt and averages them; then the calculated flux is passed to variable update; - implement subroutine advance_solution() which advances an array of conserved variables using fluxes provided as an argument; - use advance_solution() in update_solution();
2011-05-01 08:57:14 -03:00
! local variables
!
real :: dxi, dyi, dzi
Add magnetic terms to advance_solution() and update_solution().
2011-05-01 09:44:19 -03:00
#if defined MHD && defined GLM
real :: decay
#endif /* MHD & GLM */
Include forcing term in update_solution().
2011-05-01 09:24:11 -03:00
#ifdef FORCE
! local arrays
!
real, dimension( 3,im,jm,km) :: f
Add injected energy to total one in adiabatic case. - if we use forcing in the case of adiabatic equation of state, the total energy of each cell must be updated by the kinetic energy injected to this cell; otherwise the energy conservation will be violated;
2011-05-01 19:21:43 -03:00
#ifdef ADI
real, dimension( im,jm,km) :: ek, dek
#endif /* ADI */
Include forcing term in update_solution().
2011-05-01 09:24:11 -03:00
#endif /* FORCE */
Implement flux integration using RK2. - implement subroutine flux_rk2() which calculates fluxes at time t and t + dt and averages them; then the calculated flux is passed to variable update; - implement subroutine advance_solution() which advances an array of conserved variables using fluxes provided as an argument; - use advance_solution() in update_solution();
2011-05-01 08:57:14 -03:00
!
!-------------------------------------------------------------------------------
!
! prepare dxi, dyi, and dzi
!
dxi = adxi(pblock%meta%level)
dyi = adyi(pblock%meta%level)
#if NDIMS == 3
dzi = adzi(pblock%meta%level)
#endif /* NDIMS == 3 */
! perform update of conserved variables of the given block using its fluxes
!
call advance_solution(pblock%u(:,:,:,:), pblock%f(:,:,:,:,:), dxi, dyi, dzi)
Add magnetic terms to advance_solution() and update_solution().
2011-05-01 09:44:19 -03:00
#if defined MHD && defined GLM
! evolve Psi due to the source term
!
decay = exp(- alpha_p * cmax * dt / dx_min)
pblock%u(iph,:,:,:) = decay * pblock%u(iph,:,:,:)
#endif /* MHD & GLM */
Include forcing term in update_solution().
2011-05-01 09:24:11 -03:00
#ifdef FORCE
! obtain the forcing terms in real space
!
call real_forcing(pblock%meta%level, pblock%meta%xmin, pblock%meta%ymin &
, pblock%meta%zmin, f(:,:,:,:))
Add injected energy to total one in adiabatic case. - if we use forcing in the case of adiabatic equation of state, the total energy of each cell must be updated by the kinetic energy injected to this cell; otherwise the energy conservation will be violated;
2011-05-01 19:21:43 -03:00
#ifdef ADI
! calculate kinetic energy before adding the forcing term
!
ek(:,:,:) = 0.5d0 * (pblock%u(imx,:,:,:)**2 + pblock%u(imy,:,:,:)**2 &
+ pblock%u(imz,:,:,:)**2) / pblock%u(idn,:,:,:)
#endif /* ADI */
Include forcing term in update_solution().
2011-05-01 09:24:11 -03:00
! update momenta due to the forcing terms
!
pblock%u(imx,:,:,:) = pblock%u(imx,:,:,:) &
+ pblock%u(idn,:,:,:) * f(inx,:,:,:)
pblock%u(imy,:,:,:) = pblock%u(imy,:,:,:) &
+ pblock%u(idn,:,:,:) * f(iny,:,:,:)
pblock%u(imz,:,:,:) = pblock%u(imz,:,:,:) &
+ pblock%u(idn,:,:,:) * f(inz,:,:,:)
Add injected energy to total one in adiabatic case. - if we use forcing in the case of adiabatic equation of state, the total energy of each cell must be updated by the kinetic energy injected to this cell; otherwise the energy conservation will be violated;
2011-05-01 19:21:43 -03:00
#ifdef ADI
! calculate kinetic energy after adding the forcing term
!
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(:,:,:)
#endif /* ADI */
Include forcing term in update_solution().
2011-05-01 09:24:11 -03:00
#endif /* FORCE */
Implement flux integration using RK2. - implement subroutine flux_rk2() which calculates fluxes at time t and t + dt and averages them; then the calculated flux is passed to variable update; - implement subroutine advance_solution() which advances an array of conserved variables using fluxes provided as an argument; - use advance_solution() in update_solution();
2011-05-01 08:57:14 -03:00
!-------------------------------------------------------------------------------
!
end subroutine update_solution
!
!===============================================================================
!
! advance_solution: subroutine performs an one step update of the conserved
! variables array using the numerical fluxes passed as an
! argument
!
!===============================================================================
!
subroutine advance_solution(u, f, dxi, dyi, dzi)
use config , only : im, jm, km
Add magnetic terms to advance_solution() and update_solution().
2011-05-01 09:44:19 -03:00
use variables, only : nqt, nfl
use variables, only : inx, iny, inz
#ifdef MHD
use variables, only : ibx, ibz
#ifdef GLM
use scheme , only : cmax
use variables, only : iph
#endif /* GLM */
#endif /* MHD */
Implement flux integration using RK2. - implement subroutine flux_rk2() which calculates fluxes at time t and t + dt and averages them; then the calculated flux is passed to variable update; - implement subroutine advance_solution() which advances an array of conserved variables using fluxes provided as an argument; - use advance_solution() in update_solution();
2011-05-01 08:57:14 -03:00
implicit none
! input arguments
!
real, dimension( nqt,im,jm,km), intent(inout) :: u
real, dimension(NDIMS,nqt,im,jm,km), intent(in) :: f
real :: dxi, dyi, dzi
Add one step update of conserved variables. - after calculating an integrals of fluxes with a given order we performed one step updated of the conserved variables stored in data blocks using corresponding numerical flux stored in the same data block; the subroutine which performs update is called update_solution();
2011-04-30 13:14:07 -03:00
! local variables
!
integer :: i, j, k, im1, jm1, km1
real :: dhx, dhy, dhz
Add magnetic terms to advance_solution() and update_solution().
2011-05-01 09:44:19 -03:00
#if defined MHD && defined GLM
real :: ch2
#endif /* MHD & GLM */
! local arrays
!
real, dimension(nqt,im,jm,km) :: du
Add one step update of conserved variables. - after calculating an integrals of fluxes with a given order we performed one step updated of the conserved variables stored in data blocks using corresponding numerical flux stored in the same data block; the subroutine which performs update is called update_solution();
2011-04-30 13:14:07 -03:00
!
!-------------------------------------------------------------------------------
!
! prepare dxi, dyi, and dzi
!
Implement flux integration using RK2. - implement subroutine flux_rk2() which calculates fluxes at time t and t + dt and averages them; then the calculated flux is passed to variable update; - implement subroutine advance_solution() which advances an array of conserved variables using fluxes provided as an argument; - use advance_solution() in update_solution();
2011-05-01 08:57:14 -03:00
dhx = dt * dxi
dhy = dt * dyi
Add one step update of conserved variables. - after calculating an integrals of fluxes with a given order we performed one step updated of the conserved variables stored in data blocks using corresponding numerical flux stored in the same data block; the subroutine which performs update is called update_solution();
2011-04-30 13:14:07 -03:00
#if NDIMS == 3
Implement flux integration using RK2. - implement subroutine flux_rk2() which calculates fluxes at time t and t + dt and averages them; then the calculated flux is passed to variable update; - implement subroutine advance_solution() which advances an array of conserved variables using fluxes provided as an argument; - use advance_solution() in update_solution();
2011-05-01 08:57:14 -03:00
dhz = dt * dzi
Add one step update of conserved variables. - after calculating an integrals of fluxes with a given order we performed one step updated of the conserved variables stored in data blocks using corresponding numerical flux stored in the same data block; the subroutine which performs update is called update_solution();
2011-04-30 13:14:07 -03:00
#endif /* NDIMS == 3 */
Add magnetic terms to advance_solution() and update_solution().
2011-05-01 09:44:19 -03:00
! reset the increment array du
!
du(:,:,:,:) = 0.0d0
Add one step update of conserved variables. - after calculating an integrals of fluxes with a given order we performed one step updated of the conserved variables stored in data blocks using corresponding numerical flux stored in the same data block; the subroutine which performs update is called update_solution();
2011-04-30 13:14:07 -03:00
! perform update along the X direction
!
do i = 1, im
im1 = max(1, i - 1)
Add magnetic terms to advance_solution() and update_solution().
2011-05-01 09:44:19 -03:00
du(:,i,:,:) = du(:,i,:,:) + dhx * (f(inx,:,im1,:,:) - f(inx,:,i,:,:))
Add one step update of conserved variables. - after calculating an integrals of fluxes with a given order we performed one step updated of the conserved variables stored in data blocks using corresponding numerical flux stored in the same data block; the subroutine which performs update is called update_solution();
2011-04-30 13:14:07 -03:00
end do
! perform update along the Y direction
!
do j = 1, jm
jm1 = max(1, j - 1)
Add magnetic terms to advance_solution() and update_solution().
2011-05-01 09:44:19 -03:00
du(:,:,j,:) = du(:,:,j,:) + dhy * (f(iny,:,:,jm1,:) - f(iny,:,:,j,:))
Add one step update of conserved variables. - after calculating an integrals of fluxes with a given order we performed one step updated of the conserved variables stored in data blocks using corresponding numerical flux stored in the same data block; the subroutine which performs update is called update_solution();
2011-04-30 13:14:07 -03:00
end do
#if NDIMS == 3
! perform update along the Z direction
!
do k = 1, km
km1 = max(1, k - 1)
Add magnetic terms to advance_solution() and update_solution().
2011-05-01 09:44:19 -03:00
du(:,:,:,k) = du(:,:,:,k) + dhz * (f(inz,:,:,:,km1) - f(inz,:,:,:,k))
Add one step update of conserved variables. - after calculating an integrals of fluxes with a given order we performed one step updated of the conserved variables stored in data blocks using corresponding numerical flux stored in the same data block; the subroutine which performs update is called update_solution();
2011-04-30 13:14:07 -03:00
end do
#endif /* NDIMS == 3 */
Add magnetic terms to advance_solution() and update_solution().
2011-05-01 09:44:19 -03:00
! update the solution for the fluid variables
!
u( 1:nfl,:,:,:) = u( 1:nfl,:,:,:) + du( 1:nfl,:,:,:)
#ifdef MHD
! update the solution for the magnetic variables
!
u(ibx:ibz,:,:,:) = u(ibx:ibz,:,:,:) + du(ibx:ibz,:,:,:)
#ifdef GLM
! calculate c_h^2
!
ch2 = cmax * cmax
! update the solution for the scalar potential Psi
!
u(iph,:,:,:) = u(iph,:,:,:) + ch2 * du(iph,:,:,:)
#endif /* GLM */
#endif /* MHD */
Add one step update of conserved variables. - after calculating an integrals of fluxes with a given order we performed one step updated of the conserved variables stored in data blocks using corresponding numerical flux stored in the same data block; the subroutine which performs update is called update_solution();
2011-04-30 13:14:07 -03:00
!-------------------------------------------------------------------------------
!
Implement flux integration using RK2. - implement subroutine flux_rk2() which calculates fluxes at time t and t + dt and averages them; then the calculated flux is passed to variable update; - implement subroutine advance_solution() which advances an array of conserved variables using fluxes provided as an argument; - use advance_solution() in update_solution();
2011-05-01 08:57:14 -03:00
end subroutine advance_solution
Implement flux integration using Euler method. - subroutine flux_euler() integrates the numerical flux using the first order Euler method;
2011-04-30 12:53:17 -03:00
#ifdef EULER
!
!===============================================================================
!
! flux_euler: subroutine performs the first order integration of the numerical
! flux
!
!===============================================================================
!
subroutine flux_euler(pblock)
use blocks , only : block_data
use mesh , only : adxi, adyi, adzi
use scheme , only : update_flux
implicit none
! input arguments
!
type(block_data), intent(inout) :: pblock
! local variables
!
real :: dxi, dyi, dzi
!
!-------------------------------------------------------------------------------
!
! prepare dxi, dyi, and dzi
!
dxi = adxi(pblock%meta%level)
dyi = adyi(pblock%meta%level)
dzi = adzi(pblock%meta%level)
! 1st step of integration
!
call update_flux(pblock%u(:,:,:,:), pblock%f(:,:,:,:,:), dxi, dyi, dzi)
!-------------------------------------------------------------------------------
!
end subroutine flux_euler
#endif /* EULER */
Implement flux integration using RK2. - implement subroutine flux_rk2() which calculates fluxes at time t and t + dt and averages them; then the calculated flux is passed to variable update; - implement subroutine advance_solution() which advances an array of conserved variables using fluxes provided as an argument; - use advance_solution() in update_solution();
2011-05-01 08:57:14 -03:00
#ifdef RK2
!
!===============================================================================
!
! flux_rk2: subroutine performs integration of the numerical flux using
! the second order Runge-Kutta method
!
!===============================================================================
!
subroutine flux_rk2(pblock)
use blocks , only : block_data
use config , only : im, jm, km
use mesh , only : adxi, adyi, adzi
use scheme , only : update_flux
use variables, only : nqt
implicit none
! input arguments
!
type(block_data), intent(inout) :: pblock
! local variables
!
real :: dxi, dyi, dzi
! local arrays
!
real, dimension( nqt,im,jm,km) :: u
real, dimension(NDIMS,nqt,im,jm,km) :: f0, f1
!
!-------------------------------------------------------------------------------
!
! prepare dxi, dyi, and dzi
!
dxi = adxi(pblock%meta%level)
dyi = adyi(pblock%meta%level)
dzi = adzi(pblock%meta%level)
! copy the initial state to the local array u
!
u(:,:,:,:) = pblock%u(:,:,:,:)
! calculate fluxes at the moment t
!
call update_flux(u(:,:,:,:), f0(:,:,:,:,:), dxi, dyi, dzi)
! advance the solution to (t + dt) using computed fluxes in this substep
!
call advance_solution(u(:,:,:,:), f0(:,:,:,:,:), dxi, dyi, dzi)
! calculate fluxes at the moment (t + dt)
!
call update_flux(u(:,:,:,:), f1(:,:,:,:,:), dxi, dyi, dzi)
! calculate the time averaged flux
!
pblock%f(:,:,:,:,:) = 0.5d0 * (f0(:,:,:,:,:) + f1(:,:,:,:,:))
!-------------------------------------------------------------------------------
!
end subroutine flux_rk2
#endif /* RK2 */
Implement flux integration using RK3 method.
2011-05-01 10:31:27 -03:00
#ifdef RK3
!
!===============================================================================
!
! flux_rk3: subroutine performs integration of the numerical flux using
! the third order Runge-Kutta method
!
!===============================================================================
!
subroutine flux_rk3(pblock)
use blocks , only : block_data
use config , only : im, jm, km
use mesh , only : adxi, adyi, adzi
use scheme , only : update_flux
use variables, only : nqt
implicit none
! input arguments
!
type(block_data), intent(inout) :: pblock
! local variables
!
real :: dxi, dyi, dzi
! local arrays
!
real, dimension( nqt,im,jm,km) :: u
real, dimension(NDIMS,nqt,im,jm,km) :: f0, f1, f2
!
!-------------------------------------------------------------------------------
!
! prepare dxi, dyi, and dzi
!
dxi = adxi(pblock%meta%level)
dyi = adyi(pblock%meta%level)
dzi = adzi(pblock%meta%level)
! copy the initial state to the local array u
!
u(:,:,:,:) = pblock%u(:,:,:,:)
! calculate fluxes at the moment t
!
call update_flux(u(:,:,:,:), f0(:,:,:,:,:), dxi, dyi, dzi)
! advance the solution to (t + dt) using computed fluxes in this substep
!
call advance_solution(u(:,:,:,:), f0(:,:,:,:,:), dxi, dyi, dzi)
! calculate fluxes at the moment (t + dt)
!
call update_flux(u(:,:,:,:), f1(:,:,:,:,:), dxi, dyi, dzi)
! copy the initial state to the local array u
!
u(:,:,:,:) = pblock%u(:,:,:,:)
! average fluxes from t and t + dt and prepare for half step update
!
f2(:,:,:,:,:) = 0.25d0 * (f0(:,:,:,:,:) + f1(:,:,:,:,:))
! advance the solution to (t + dt/2) using computed flux
!
call advance_solution(u(:,:,:,:), f2(:,:,:,:,:), dxi, dyi, dzi)
! calculate fluxes at the moment (t + dt/2)
!
call update_flux(u(:,:,:,:), f2(:,:,:,:,:), dxi, dyi, dzi)
! calculate the time averaged flux using Gauss formula
!
pblock%f(:,:,:,:,:) = (f0(:,:,:,:,:) + f1(:,:,:,:,:) &
+ 4.0d0 * f2(:,:,:,:,:)) / 6.0d0
!-------------------------------------------------------------------------------
!
end subroutine flux_rk3
#endif /* RK3 */
Add a subroutine to calculate numerical fluxes. - the subroutine update_flux() in the SCHEME module calculates numerical fluxes from the conserved variables; the calculated fluxes are then returned;
2011-04-30 12:28:02 -03:00
#else /* CONSERVATIVE */
New time integration method EULER.
2008-12-09 21:06:21 -06:00
#ifdef EULER
!
!===============================================================================
!
! evolve_euler: subroutine evolves the current block using Euler integration
!
!===============================================================================
!
subroutine evolve_euler(pblock)
Move variable indices to new module 'variables'. VARIABLES - create new module 'variables' which stores references to variable indices; we gonna store dofferent objects related to variables in this module;
2010-12-01 09:25:30 -02:00
use blocks , only : block_data
use config , only : im, jm, km
EVOLUTION: add forcing terms in the block update.
2010-12-08 21:48:43 -02:00
#ifdef FORCE
use forcing , only : real_forcing
#endif /* FORCE */
Revert "INTERPOLATION: pass the spacial increment to reconstruct()." This reverts commit d92ae5d3a62743b66e97c42be4aa5f31ef5993d1.
2010-12-06 16:20:49 -02:00
use mesh , only : adxi, adyi, adzi
Add support for the shapes in the domain. Now, we can define any shape for a given problem inside the domain, which is not updated during the evolution. This allows for using the sources of any kind in the problem studies, such as the colliding winds in the binary stars.
2008-12-18 12:18:36 -06:00
#ifdef SHAPE
Move variable indices to new module 'variables'. VARIABLES - create new module 'variables' which stores references to variable indices; we gonna store dofferent objects related to variables in this module;
2010-12-01 09:25:30 -02:00
use problem , only : update_shapes
Add support for the shapes in the domain. Now, we can define any shape for a given problem inside the domain, which is not updated during the evolution. This allows for using the sources of any kind in the problem studies, such as the colliding winds in the binary stars.
2008-12-18 12:18:36 -06:00
#endif /* SHAPE */
SCHEME: move cmax parameter from module 'evolution' to 'scheme'.
2010-12-01 13:13:27 -02:00
use scheme , only : update, cmax
Move variable indices to new module 'variables'. VARIABLES - create new module 'variables' which stores references to variable indices; we gonna store dofferent objects related to variables in this module;
2010-12-01 09:25:30 -02:00
use variables, only : nqt, nfl
GLM-MHD: update global solution. - update the magnetic field components and scalar potential in the Euler and RK2 methods, evolve_euler() and evolve_rk2() subroutines, respectively;
2010-12-01 10:53:21 -02:00
#ifdef MHD
use variables, only : ibx, ibz
#ifdef GLM
GLM-MHD: implement Psi decay due to a source term. - the scalar potential Psi evolution is controlled by a dissipative source term; include this update using a simple analytical solution;
2010-12-01 11:20:25 -02:00
use config , only : alpha_p
use mesh , only : dx_min
GLM-MHD: update global solution. - update the magnetic field components and scalar potential in the Euler and RK2 methods, evolve_euler() and evolve_rk2() subroutines, respectively;
2010-12-01 10:53:21 -02:00
use variables, only : iph
#endif /* GLM */
#endif /* MHD */
EVOLUTION: add forcing terms in the block update.
2010-12-08 21:48:43 -02:00
#ifdef FORCE
use variables, only : idn, imx, imy, imz
Add injected energy to total one in adiabatic case. - if we use forcing in the case of adiabatic equation of state, the total energy of each cell must be updated by the kinetic energy injected to this cell; otherwise the energy conservation will be violated;
2011-05-01 19:21:43 -03:00
#ifdef ADI
use variables, only : ien
#endif /* ADI */
EVOLUTION: add forcing terms in the block update.
2010-12-08 21:48:43 -02:00
#endif /* FORCE */
New time integration method EULER.
2008-12-09 21:06:21 -06:00
implicit none
! input arguments
!
Initial support for the MHD equations. VARIABLES - add indices for the magnetic field components, both face and cell centered SOLVER, MHD - add support for the magnetohydrodynamic (MHD) equations to the subroutines cons2prim(), prim2cons() - add MHD flux and the fastest speed calculation in the subroutine fluxspeed() - include magnetosonic speed in the calculation of the maximum speed in the system required for estimation of the new time step - extend the HLL solver in subroutine hll() to support MHD - calculate the magnetic field update according to a CT scheme in the subroutine update() INTERPOLATION - add subroutine magtocen() to interpolate cell centered magnetic field EVOLUTION - add evolution of the magnetic field components in the evolve_euler() and evolve_rk2() time integration subroutines - also call the subroutine magtocen() in the right places BOUNDARY CONDITIONS - support for magnetic field boundary update only in the case of blocks with the same level so far; later we need to include proper restriction and prolongation for the magnetic field to keep its divergence equal zero PROBLEMS, BLAST - extend the blast problem to include the initial magnetic field IO, HDF5 - write down in a HDF5 file magnetic field components
2009-10-28 00:12:18 -02:00
type(block_data), intent(inout) :: pblock
New time integration method EULER.
2008-12-09 21:06:21 -06:00
! local variables
!
Revert "INTERPOLATION: pass the spacial increment to reconstruct()." This reverts commit d92ae5d3a62743b66e97c42be4aa5f31ef5993d1.
2010-12-06 16:20:49 -02:00
real :: dxi, dyi, dzi, ch2
GLM-MHD: implement Psi decay due to a source term. - the scalar potential Psi evolution is controlled by a dissipative source term; include this update using a simple analytical solution;
2010-12-01 11:20:25 -02:00
#if defined MHD && defined GLM
Revert "INTERPOLATION: pass the spacial increment to reconstruct()." This reverts commit d92ae5d3a62743b66e97c42be4aa5f31ef5993d1.
2010-12-06 16:20:49 -02:00
real :: decay
GLM-MHD: implement Psi decay due to a source term. - the scalar potential Psi evolution is controlled by a dissipative source term; include this update using a simple analytical solution;
2010-12-01 11:20:25 -02:00
#endif /* MHD & GLM */
New time integration method EULER.
2008-12-09 21:06:21 -06:00
! local arrays
!
Remove flux calculation from the subroutine update(). SCHEME - remove the flux calculation from the subroutine update() and all dependent subroutines;
2010-09-19 00:08:20 +02:00
real, dimension(nqt,im,jm,km) :: du
EVOLUTION: add forcing terms in the block update.
2010-12-08 21:48:43 -02:00
#ifdef FORCE
real, dimension( 3,im,jm,km) :: f
Add injected energy to total one in adiabatic case. - if we use forcing in the case of adiabatic equation of state, the total energy of each cell must be updated by the kinetic energy injected to this cell; otherwise the energy conservation will be violated;
2011-05-01 19:21:43 -03:00
#ifdef ADI
real, dimension( im,jm,km) :: ek, dek
#endif /* ADI */
EVOLUTION: add forcing terms in the block update.
2010-12-08 21:48:43 -02:00
#endif /* FORCE */
New time integration method EULER.
2008-12-09 21:06:21 -06:00
!
!-------------------------------------------------------------------------------
!
Revert "INTERPOLATION: pass the spacial increment to reconstruct()." This reverts commit d92ae5d3a62743b66e97c42be4aa5f31ef5993d1.
2010-12-06 16:20:49 -02:00
! prepare dxi, dyi, and dzi
!
dxi = adxi(pblock%meta%level)
dyi = adyi(pblock%meta%level)
dzi = adzi(pblock%meta%level)
FORCE: update momenta due to forcinf only once. - do not put forcing terms in the RK substeps; instead, update the momenta components only once before executing the Riemann solver and update step;
2011-03-06 23:37:51 -03:00
! 1st step of integration
!
call update(pblock%u(:,:,:,:), du(:,:,:,:), dxi, dyi, dzi)
Add support for the shapes in the domain. Now, we can define any shape for a given problem inside the domain, which is not updated during the evolution. This allows for using the sources of any kind in the problem studies, such as the colliding winds in the binary stars.
2008-12-18 12:18:36 -06:00
#ifdef SHAPE
! restrict update in a defined shape
!
Revert "INTERPOLATION: pass the spacial increment to reconstruct()." This reverts commit d92ae5d3a62743b66e97c42be4aa5f31ef5993d1.
2010-12-06 16:20:49 -02:00
call update_shapes(pblock, du)
Add support for the shapes in the domain. Now, we can define any shape for a given problem inside the domain, which is not updated during the evolution. This allows for using the sources of any kind in the problem studies, such as the colliding winds in the binary stars.
2008-12-18 12:18:36 -06:00
#endif /* SHAPE */
New time integration method EULER.
2008-12-09 21:06:21 -06:00
GLM-MHD: update global solution. - update the magnetic field components and scalar potential in the Euler and RK2 methods, evolve_euler() and evolve_rk2() subroutines, respectively;
2010-12-01 10:53:21 -02:00
! update the solution for the fluid variables
New time integration method EULER.
2008-12-09 21:06:21 -06:00
!
GLM-MHD: update global solution. - update the magnetic field components and scalar potential in the Euler and RK2 methods, evolve_euler() and evolve_rk2() subroutines, respectively;
2010-12-01 10:53:21 -02:00
pblock%u(1:nfl,:,:,:) = pblock%u(1:nfl,:,:,:) + dt * du(1:nfl,:,:,:)
#ifdef MHD
! update the solution for the magnetic variables
!
pblock%u(ibx:ibz,:,:,:) = pblock%u(ibx:ibz,:,:,:) + dt * du(ibx:ibz,:,:,:)
#ifdef GLM
GLM-MHD: take into account the maximum speed in update. - the divergence of B propagates with the maximum speed c_h; take it into account while updating the solution for the scalar potential Psi;
2010-12-01 10:57:40 -02:00
! calculate c_h^2
!
ch2 = cmax * cmax
GLM-MHD: update global solution. - update the magnetic field components and scalar potential in the Euler and RK2 methods, evolve_euler() and evolve_rk2() subroutines, respectively;
2010-12-01 10:53:21 -02:00
! update the solution for the scalar potential Psi
!
GLM-MHD: take into account the maximum speed in update. - the divergence of B propagates with the maximum speed c_h; take it into account while updating the solution for the scalar potential Psi;
2010-12-01 10:57:40 -02:00
pblock%u(iph,:,:,:) = pblock%u(iph,:,:,:) + ch2 * dt * du(iph,:,:,:)
GLM-MHD: implement Psi decay due to a source term. - the scalar potential Psi evolution is controlled by a dissipative source term; include this update using a simple analytical solution;
2010-12-01 11:20:25 -02:00
! evolve Psi due to the source term
!
decay = exp(- alpha_p * cmax * dt / dx_min)
pblock%u(iph,:,:,:) = decay * pblock%u(iph,:,:,:)
GLM-MHD: update global solution. - update the magnetic field components and scalar potential in the Euler and RK2 methods, evolve_euler() and evolve_rk2() subroutines, respectively;
2010-12-01 10:53:21 -02:00
#endif /* GLM */
#endif /* MHD */
Add forcing terms after the main update. - change order of calculating terms; first perform the main update of the conserved variables, then add forcing and source terms; in this way we do not change time step which determines the stability of scheme before the conserved variables enter Riemann solver; this is especially important in highly supersonic simulations, when in the presence of very low density, even small change of velocity can make the scheme numerically unstable;
2011-04-26 10:36:16 -03:00
#ifdef FORCE
! obtain the forcing terms in real space
!
call real_forcing(pblock%meta%level, pblock%meta%xmin, pblock%meta%ymin &
, pblock%meta%zmin, f(:,:,:,:))
Add injected energy to total one in adiabatic case. - if we use forcing in the case of adiabatic equation of state, the total energy of each cell must be updated by the kinetic energy injected to this cell; otherwise the energy conservation will be violated;
2011-05-01 19:21:43 -03:00
#ifdef ADI
! calculate kinetic energy before adding the forcing term
!
ek(:,:,:) = 0.5d0 * (pblock%u(imx,:,:,:)**2 + pblock%u(imy,:,:,:)**2 &
+ pblock%u(imz,:,:,:)**2) / pblock%u(idn,:,:,:)
#endif /* ADI */
Add forcing terms after the main update. - change order of calculating terms; first perform the main update of the conserved variables, then add forcing and source terms; in this way we do not change time step which determines the stability of scheme before the conserved variables enter Riemann solver; this is especially important in highly supersonic simulations, when in the presence of very low density, even small change of velocity can make the scheme numerically unstable;
2011-04-26 10:36:16 -03:00
! 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,:,:,:)
Add injected energy to total one in adiabatic case. - if we use forcing in the case of adiabatic equation of state, the total energy of each cell must be updated by the kinetic energy injected to this cell; otherwise the energy conservation will be violated;
2011-05-01 19:21:43 -03:00
#ifdef ADI
! calculate kinetic energy after adding the forcing term
!
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(:,:,:)
#endif /* ADI */
Add forcing terms after the main update. - change order of calculating terms; first perform the main update of the conserved variables, then add forcing and source terms; in this way we do not change time step which determines the stability of scheme before the conserved variables enter Riemann solver; this is especially important in highly supersonic simulations, when in the presence of very low density, even small change of velocity can make the scheme numerically unstable;
2011-04-26 10:36:16 -03:00
#endif /* FORCE */
Initial support for the MHD equations. VARIABLES - add indices for the magnetic field components, both face and cell centered SOLVER, MHD - add support for the magnetohydrodynamic (MHD) equations to the subroutines cons2prim(), prim2cons() - add MHD flux and the fastest speed calculation in the subroutine fluxspeed() - include magnetosonic speed in the calculation of the maximum speed in the system required for estimation of the new time step - extend the HLL solver in subroutine hll() to support MHD - calculate the magnetic field update according to a CT scheme in the subroutine update() INTERPOLATION - add subroutine magtocen() to interpolate cell centered magnetic field EVOLUTION - add evolution of the magnetic field components in the evolve_euler() and evolve_rk2() time integration subroutines - also call the subroutine magtocen() in the right places BOUNDARY CONDITIONS - support for magnetic field boundary update only in the case of blocks with the same level so far; later we need to include proper restriction and prolongation for the magnetic field to keep its divergence equal zero PROBLEMS, BLAST - extend the blast problem to include the initial magnetic field IO, HDF5 - write down in a HDF5 file magnetic field components
2009-10-28 00:12:18 -02:00
!
New time integration method EULER.
2008-12-09 21:06:21 -06:00
!-------------------------------------------------------------------------------
!
end subroutine evolve_euler
#endif /* EULER */
Runge-Kutta 2nd order time integration implemented. In addition, I've done some fixes to the problem initialization, and I defined new variables igrids, jgrids, kgrids, which specify the dimensions of the block.
2008-12-08 16:21:59 -06:00
#ifdef RK2
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
!
Updated comments in the evolution module.
2008-12-08 21:07:10 -06:00
!===============================================================================
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
!
EVOLUTION: implement the 3rd order RK method. - implement the 3rd order Runge-Kutta time integration; - rewrite slightly the 2nd order Runge-Kutta method;
2010-12-03 15:38:48 -02:00
! evolve_rk2: subroutine evolves the current block using the 2nd order
! Runge-Kutta method
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
!
Updated comments in the evolution module.
2008-12-08 21:07:10 -06:00
!===============================================================================
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
!
subroutine evolve_rk2(pblock)
Move variable indices to new module 'variables'. VARIABLES - create new module 'variables' which stores references to variable indices; we gonna store dofferent objects related to variables in this module;
2010-12-01 09:25:30 -02:00
use blocks , only : block_data
use config , only : im, jm, km
EVOLUTION: add forcing terms in the RK2 update.
2010-12-08 21:55:45 -02:00
#ifdef FORCE
use forcing , only : real_forcing
#endif /* FORCE */
Revert "INTERPOLATION: pass the spacial increment to reconstruct()." This reverts commit d92ae5d3a62743b66e97c42be4aa5f31ef5993d1.
2010-12-06 16:20:49 -02:00
use mesh , only : adxi, adyi, adzi
Add support for the shapes in the domain. Now, we can define any shape for a given problem inside the domain, which is not updated during the evolution. This allows for using the sources of any kind in the problem studies, such as the colliding winds in the binary stars.
2008-12-18 12:18:36 -06:00
#ifdef SHAPE
Move variable indices to new module 'variables'. VARIABLES - create new module 'variables' which stores references to variable indices; we gonna store dofferent objects related to variables in this module;
2010-12-01 09:25:30 -02:00
use problem , only : update_shapes
Add support for the shapes in the domain. Now, we can define any shape for a given problem inside the domain, which is not updated during the evolution. This allows for using the sources of any kind in the problem studies, such as the colliding winds in the binary stars.
2008-12-18 12:18:36 -06:00
#endif /* SHAPE */
SCHEME: move cmax parameter from module 'evolution' to 'scheme'.
2010-12-01 13:13:27 -02:00
use scheme , only : update, cmax
Move variable indices to new module 'variables'. VARIABLES - create new module 'variables' which stores references to variable indices; we gonna store dofferent objects related to variables in this module;
2010-12-01 09:25:30 -02:00
use variables, only : nqt, nfl
GLM-MHD: update global solution. - update the magnetic field components and scalar potential in the Euler and RK2 methods, evolve_euler() and evolve_rk2() subroutines, respectively;
2010-12-01 10:53:21 -02:00
#ifdef MHD
use variables, only : ibx, ibz
#ifdef GLM
GLM-MHD: implement Psi decay due to a source term. - the scalar potential Psi evolution is controlled by a dissipative source term; include this update using a simple analytical solution;
2010-12-01 11:20:25 -02:00
use config , only : alpha_p
use mesh , only : dx_min
GLM-MHD: update global solution. - update the magnetic field components and scalar potential in the Euler and RK2 methods, evolve_euler() and evolve_rk2() subroutines, respectively;
2010-12-01 10:53:21 -02:00
use variables, only : iph
#endif /* GLM */
#endif /* MHD */
EVOLUTION: add forcing terms in the RK2 update.
2010-12-08 21:55:45 -02:00
#ifdef FORCE
use variables, only : idn, imx, imy, imz
Add injected energy to total one in adiabatic case. - if we use forcing in the case of adiabatic equation of state, the total energy of each cell must be updated by the kinetic energy injected to this cell; otherwise the energy conservation will be violated;
2011-05-01 19:21:43 -03:00
#ifdef ADI
use variables, only : ien
#endif /* ADI */
EVOLUTION: add forcing terms in the RK2 update.
2010-12-08 21:55:45 -02:00
#endif /* FORCE */
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
implicit none
! input arguments
!
Initial support for the MHD equations. VARIABLES - add indices for the magnetic field components, both face and cell centered SOLVER, MHD - add support for the magnetohydrodynamic (MHD) equations to the subroutines cons2prim(), prim2cons() - add MHD flux and the fastest speed calculation in the subroutine fluxspeed() - include magnetosonic speed in the calculation of the maximum speed in the system required for estimation of the new time step - extend the HLL solver in subroutine hll() to support MHD - calculate the magnetic field update according to a CT scheme in the subroutine update() INTERPOLATION - add subroutine magtocen() to interpolate cell centered magnetic field EVOLUTION - add evolution of the magnetic field components in the evolve_euler() and evolve_rk2() time integration subroutines - also call the subroutine magtocen() in the right places BOUNDARY CONDITIONS - support for magnetic field boundary update only in the case of blocks with the same level so far; later we need to include proper restriction and prolongation for the magnetic field to keep its divergence equal zero PROBLEMS, BLAST - extend the blast problem to include the initial magnetic field IO, HDF5 - write down in a HDF5 file magnetic field components
2009-10-28 00:12:18 -02:00
type(block_data), intent(inout) :: pblock
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
! local variables
!
Revert "INTERPOLATION: pass the spacial increment to reconstruct()." This reverts commit d92ae5d3a62743b66e97c42be4aa5f31ef5993d1.
2010-12-06 16:20:49 -02:00
real :: dxi, dyi, dzi, ch2
GLM-MHD: implement Psi decay due to a source term. - the scalar potential Psi evolution is controlled by a dissipative source term; include this update using a simple analytical solution;
2010-12-01 11:20:25 -02:00
#if defined MHD && defined GLM
Revert "INTERPOLATION: pass the spacial increment to reconstruct()." This reverts commit d92ae5d3a62743b66e97c42be4aa5f31ef5993d1.
2010-12-06 16:20:49 -02:00
real :: decay
GLM-MHD: implement Psi decay due to a source term. - the scalar potential Psi evolution is controlled by a dissipative source term; include this update using a simple analytical solution;
2010-12-01 11:20:25 -02:00
#endif /* MHD & GLM */
Runge-Kutta 2nd order time integration implemented. In addition, I've done some fixes to the problem initialization, and I defined new variables igrids, jgrids, kgrids, which specify the dimensions of the block.
2008-12-08 16:21:59 -06:00
! local arrays
!
Remove flux calculation from the subroutine update(). SCHEME - remove the flux calculation from the subroutine update() and all dependent subroutines;
2010-09-19 00:08:20 +02:00
real, dimension(nqt,im,jm,km) :: u1, du
EVOLUTION: add forcing terms in the RK2 update.
2010-12-08 21:55:45 -02:00
#ifdef FORCE
real, dimension( 3,im,jm,km) :: f
Add injected energy to total one in adiabatic case. - if we use forcing in the case of adiabatic equation of state, the total energy of each cell must be updated by the kinetic energy injected to this cell; otherwise the energy conservation will be violated;
2011-05-01 19:21:43 -03:00
#ifdef ADI
real, dimension( im,jm,km) :: ek, dek
#endif /* ADI */
EVOLUTION: add forcing terms in the RK2 update.
2010-12-08 21:55:45 -02:00
#endif /* FORCE */
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
!
Updated comments in the evolution module.
2008-12-08 21:07:10 -06:00
!-------------------------------------------------------------------------------
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
!
Revert "INTERPOLATION: pass the spacial increment to reconstruct()." This reverts commit d92ae5d3a62743b66e97c42be4aa5f31ef5993d1.
2010-12-06 16:20:49 -02:00
! prepare dxi, dyi, and dzi
!
dxi = adxi(pblock%meta%level)
dyi = adyi(pblock%meta%level)
dzi = adzi(pblock%meta%level)
EVOLUTION: implement the 3rd order RK method. - implement the 3rd order Runge-Kutta time integration; - rewrite slightly the 2nd order Runge-Kutta method;
2010-12-03 15:38:48 -02:00
!! 1st step of integration
!!
Revert "INTERPOLATION: pass the spacial increment to reconstruct()." This reverts commit d92ae5d3a62743b66e97c42be4aa5f31ef5993d1.
2010-12-06 16:20:49 -02:00
call update(pblock%u(:,:,:,:), du(:,:,:,:), dxi, dyi, dzi)
Add support for the shapes in the domain. Now, we can define any shape for a given problem inside the domain, which is not updated during the evolution. This allows for using the sources of any kind in the problem studies, such as the colliding winds in the binary stars.
2008-12-18 12:18:36 -06:00
#ifdef SHAPE
! restrict update in a defined shape
!
EVOLUTION: implement the 3rd order RK method. - implement the 3rd order Runge-Kutta time integration; - rewrite slightly the 2nd order Runge-Kutta method;
2010-12-03 15:38:48 -02:00
call update_shapes(pblock, du(:,:,:,:))
Runge-Kutta 2nd order time integration implemented. In addition, I've done some fixes to the problem initialization, and I defined new variables igrids, jgrids, kgrids, which specify the dimensions of the block.
2008-12-08 16:21:59 -06:00
Fix forcing term in subroutine evolve_rk2().
2011-04-26 11:50:57 -03:00
#endif /* SHAPE */
GLM-MHD: update global solution. - update the magnetic field components and scalar potential in the Euler and RK2 methods, evolve_euler() and evolve_rk2() subroutines, respectively;
2010-12-01 10:53:21 -02:00
! update the solution for the fluid variables
!
u1(1:nfl,:,:,:) = pblock%u(1:nfl,:,:,:) + dt * du(1:nfl,:,:,:)
#ifdef MHD
! update the solution for the magnetic variables
Runge-Kutta 2nd order time integration implemented. In addition, I've done some fixes to the problem initialization, and I defined new variables igrids, jgrids, kgrids, which specify the dimensions of the block.
2008-12-08 16:21:59 -06:00
!
GLM-MHD: update global solution. - update the magnetic field components and scalar potential in the Euler and RK2 methods, evolve_euler() and evolve_rk2() subroutines, respectively;
2010-12-01 10:53:21 -02:00
u1(ibx:ibz,:,:,:) = pblock%u(ibx:ibz,:,:,:) + dt * du(ibx:ibz,:,:,:)
#ifdef GLM
FORCE: update momenta due to forcinf only once. - do not put forcing terms in the RK substeps; instead, update the momenta components only once before executing the Riemann solver and update step;
2011-03-06 23:37:51 -03:00
! calculate c_h^2
!
ch2 = cmax * cmax
GLM-MHD: update global solution. - update the magnetic field components and scalar potential in the Euler and RK2 methods, evolve_euler() and evolve_rk2() subroutines, respectively;
2010-12-01 10:53:21 -02:00
! update the solution for the scalar potential Psi
!
GLM-MHD: take into account the maximum speed in update. - the divergence of B propagates with the maximum speed c_h; take it into account while updating the solution for the scalar potential Psi;
2010-12-01 10:57:40 -02:00
u1(iph,:,:,:) = pblock%u(iph,:,:,:) + ch2 * dt * du(iph,:,:,:)
Fix forcing term in subroutine evolve_rk2().
2011-04-26 11:50:57 -03:00
GLM-MHD: update global solution. - update the magnetic field components and scalar potential in the Euler and RK2 methods, evolve_euler() and evolve_rk2() subroutines, respectively;
2010-12-01 10:53:21 -02:00
#endif /* GLM */
#endif /* MHD */
Runge-Kutta 2nd order time integration implemented. In addition, I've done some fixes to the problem initialization, and I defined new variables igrids, jgrids, kgrids, which specify the dimensions of the block.
2008-12-08 16:21:59 -06:00
! 2nd step of integration
!
Revert "INTERPOLATION: pass the spacial increment to reconstruct()." This reverts commit d92ae5d3a62743b66e97c42be4aa5f31ef5993d1.
2010-12-06 16:20:49 -02:00
call update(u1(:,:,:,:), du(:,:,:,:), dxi, dyi, dzi)
Runge-Kutta 2nd order time integration implemented. In addition, I've done some fixes to the problem initialization, and I defined new variables igrids, jgrids, kgrids, which specify the dimensions of the block.
2008-12-08 16:21:59 -06:00
Add support for the shapes in the domain. Now, we can define any shape for a given problem inside the domain, which is not updated during the evolution. This allows for using the sources of any kind in the problem studies, such as the colliding winds in the binary stars.
2008-12-18 12:18:36 -06:00
#ifdef SHAPE
! restrict update in a defined shape
!
EVOLUTION: implement the 3rd order RK method. - implement the 3rd order Runge-Kutta time integration; - rewrite slightly the 2nd order Runge-Kutta method;
2010-12-03 15:38:48 -02:00
call update_shapes(pblock, du(:,:,:,:))
Add support for the shapes in the domain. Now, we can define any shape for a given problem inside the domain, which is not updated during the evolution. This allows for using the sources of any kind in the problem studies, such as the colliding winds in the binary stars.
2008-12-18 12:18:36 -06:00
Fix forcing term in subroutine evolve_rk2().
2011-04-26 11:50:57 -03:00
#endif /* SHAPE */
GLM-MHD: update global solution. - update the magnetic field components and scalar potential in the Euler and RK2 methods, evolve_euler() and evolve_rk2() subroutines, respectively;
2010-12-01 10:53:21 -02:00
! update the solution for the fluid variables
!
pblock%u(1:nfl,:,:,:) = 0.5d0 * (pblock%u(1:nfl,:,:,:) &
EVOLUTION: implement the 3rd order RK method. - implement the 3rd order Runge-Kutta time integration; - rewrite slightly the 2nd order Runge-Kutta method;
2010-12-03 15:38:48 -02:00
+ u1(1:nfl,:,:,:) + dt * du(1:nfl,:,:,:))
GLM-MHD: update global solution. - update the magnetic field components and scalar potential in the Euler and RK2 methods, evolve_euler() and evolve_rk2() subroutines, respectively;
2010-12-01 10:53:21 -02:00
#ifdef MHD
! update the solution for the magnetic variables
!
pblock%u(ibx:ibz,:,:,:) = 0.5d0 * (pblock%u(ibx:ibz,:,:,:) &
EVOLUTION: implement the 3rd order RK method. - implement the 3rd order Runge-Kutta time integration; - rewrite slightly the 2nd order Runge-Kutta method;
2010-12-03 15:38:48 -02:00
+ u1(ibx:ibz,:,:,:) + dt * du(ibx:ibz,:,:,:))
GLM-MHD: update global solution. - update the magnetic field components and scalar potential in the Euler and RK2 methods, evolve_euler() and evolve_rk2() subroutines, respectively;
2010-12-01 10:53:21 -02:00
#ifdef GLM
! update the solution for the scalar potential Psi
Runge-Kutta 2nd order time integration implemented. In addition, I've done some fixes to the problem initialization, and I defined new variables igrids, jgrids, kgrids, which specify the dimensions of the block.
2008-12-08 16:21:59 -06:00
!
GLM-MHD: update global solution. - update the magnetic field components and scalar potential in the Euler and RK2 methods, evolve_euler() and evolve_rk2() subroutines, respectively;
2010-12-01 10:53:21 -02:00
pblock%u(iph,:,:,:) = 0.5d0 * (pblock%u(iph,:,:,:) &
EVOLUTION: implement the 3rd order RK method. - implement the 3rd order Runge-Kutta time integration; - rewrite slightly the 2nd order Runge-Kutta method;
2010-12-03 15:38:48 -02:00
+ u1(iph,:,:,:) + ch2 * dt * du(iph,:,:,:))
GLM-MHD: implement Psi decay due to a source term. - the scalar potential Psi evolution is controlled by a dissipative source term; include this update using a simple analytical solution;
2010-12-01 11:20:25 -02:00
! evolve Psi due to the source term
!
decay = exp(- alpha_p * cmax * dt / dx_min)
pblock%u(iph,:,:,:) = decay * pblock%u(iph,:,:,:)
Fix forcing term in subroutine evolve_rk2().
2011-04-26 11:50:57 -03:00
GLM-MHD: update global solution. - update the magnetic field components and scalar potential in the Euler and RK2 methods, evolve_euler() and evolve_rk2() subroutines, respectively;
2010-12-01 10:53:21 -02:00
#endif /* GLM */
#endif /* MHD */
Add forcing terms after the main update. - change order of calculating terms; first perform the main update of the conserved variables, then add forcing and source terms; in this way we do not change time step which determines the stability of scheme before the conserved variables enter Riemann solver; this is especially important in highly supersonic simulations, when in the presence of very low density, even small change of velocity can make the scheme numerically unstable;
2011-04-26 10:36:16 -03:00
#ifdef FORCE
! obtain the forcing terms in real space
!
call real_forcing(pblock%meta%level, pblock%meta%xmin, pblock%meta%ymin &
, pblock%meta%zmin, f(:,:,:,:))
Add injected energy to total one in adiabatic case. - if we use forcing in the case of adiabatic equation of state, the total energy of each cell must be updated by the kinetic energy injected to this cell; otherwise the energy conservation will be violated;
2011-05-01 19:21:43 -03:00
#ifdef ADI
! calculate kinetic energy before adding the forcing term
!
ek(:,:,:) = 0.5d0 * (pblock%u(imx,:,:,:)**2 + pblock%u(imy,:,:,:)**2 &
+ pblock%u(imz,:,:,:)**2) / pblock%u(idn,:,:,:)
#endif /* ADI */
Add forcing terms after the main update. - change order of calculating terms; first perform the main update of the conserved variables, then add forcing and source terms; in this way we do not change time step which determines the stability of scheme before the conserved variables enter Riemann solver; this is especially important in highly supersonic simulations, when in the presence of very low density, even small change of velocity can make the scheme numerically unstable;
2011-04-26 10:36:16 -03:00
! 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,:,:,:)
Add injected energy to total one in adiabatic case. - if we use forcing in the case of adiabatic equation of state, the total energy of each cell must be updated by the kinetic energy injected to this cell; otherwise the energy conservation will be violated;
2011-05-01 19:21:43 -03:00
#ifdef ADI
! calculate kinetic energy after adding the forcing term
!
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(:,:,:)
#endif /* ADI */
Add forcing terms after the main update. - change order of calculating terms; first perform the main update of the conserved variables, then add forcing and source terms; in this way we do not change time step which determines the stability of scheme before the conserved variables enter Riemann solver; this is especially important in highly supersonic simulations, when in the presence of very low density, even small change of velocity can make the scheme numerically unstable;
2011-04-26 10:36:16 -03:00
#endif /* FORCE */
Initial support for the MHD equations. VARIABLES - add indices for the magnetic field components, both face and cell centered SOLVER, MHD - add support for the magnetohydrodynamic (MHD) equations to the subroutines cons2prim(), prim2cons() - add MHD flux and the fastest speed calculation in the subroutine fluxspeed() - include magnetosonic speed in the calculation of the maximum speed in the system required for estimation of the new time step - extend the HLL solver in subroutine hll() to support MHD - calculate the magnetic field update according to a CT scheme in the subroutine update() INTERPOLATION - add subroutine magtocen() to interpolate cell centered magnetic field EVOLUTION - add evolution of the magnetic field components in the evolve_euler() and evolve_rk2() time integration subroutines - also call the subroutine magtocen() in the right places BOUNDARY CONDITIONS - support for magnetic field boundary update only in the case of blocks with the same level so far; later we need to include proper restriction and prolongation for the magnetic field to keep its divergence equal zero PROBLEMS, BLAST - extend the blast problem to include the initial magnetic field IO, HDF5 - write down in a HDF5 file magnetic field components
2009-10-28 00:12:18 -02:00
!
Updated comments in the evolution module.
2008-12-08 21:07:10 -06:00
!-------------------------------------------------------------------------------
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
!
end subroutine evolve_rk2
Runge-Kutta 2nd order time integration implemented. In addition, I've done some fixes to the problem initialization, and I defined new variables igrids, jgrids, kgrids, which specify the dimensions of the block.
2008-12-08 16:21:59 -06:00
#endif /* RK2 */
EVOLUTION: implement the 3rd order RK method. - implement the 3rd order Runge-Kutta time integration; - rewrite slightly the 2nd order Runge-Kutta method;
2010-12-03 15:38:48 -02:00
#ifdef RK3
!
!===============================================================================
!
! evolve_rk3: subroutine evolves the current block using the 3rd order
! Runge-Kutta method
!
!===============================================================================
!
subroutine evolve_rk3(pblock)
use blocks , only : block_data
use config , only : im, jm, km
EVOLUTION: include forcing terms in the RK3 integration.
2010-12-09 17:54:59 -02:00
#ifdef FORCE
use forcing , only : real_forcing
#endif /* FORCE */
Revert "INTERPOLATION: pass the spacial increment to reconstruct()." This reverts commit d92ae5d3a62743b66e97c42be4aa5f31ef5993d1.
2010-12-06 16:20:49 -02:00
use mesh , only : adxi, adyi, adzi
EVOLUTION: implement the 3rd order RK method. - implement the 3rd order Runge-Kutta time integration; - rewrite slightly the 2nd order Runge-Kutta method;
2010-12-03 15:38:48 -02:00
#ifdef SHAPE
use problem , only : update_shapes
#endif /* SHAPE */
use scheme , only : update, cmax
use variables, only : nqt, nfl
#ifdef MHD
use variables, only : ibx, ibz
#ifdef GLM
use config , only : alpha_p
use mesh , only : dx_min
use variables, only : iph
#endif /* GLM */
#endif /* MHD */
EVOLUTION: include forcing terms in the RK3 integration.
2010-12-09 17:54:59 -02:00
#ifdef FORCE
use variables, only : idn, imx, imy, imz
Add injected energy to total one in adiabatic case. - if we use forcing in the case of adiabatic equation of state, the total energy of each cell must be updated by the kinetic energy injected to this cell; otherwise the energy conservation will be violated;
2011-05-01 19:21:43 -03:00
#ifdef ADI
use variables, only : ien
#endif /* ADI */
EVOLUTION: include forcing terms in the RK3 integration.
2010-12-09 17:54:59 -02:00
#endif /* FORCE */
EVOLUTION: implement the 3rd order RK method. - implement the 3rd order Runge-Kutta time integration; - rewrite slightly the 2nd order Runge-Kutta method;
2010-12-03 15:38:48 -02:00
implicit none
! input arguments
!
type(block_data), intent(inout) :: pblock
! local variables
!
Add magnetic terms to advance_solution() and update_solution().
2011-05-01 09:44:19 -03:00
real :: dxi, dyi, dzi
EVOLUTION: implement the 3rd order RK method. - implement the 3rd order Runge-Kutta time integration; - rewrite slightly the 2nd order Runge-Kutta method;
2010-12-03 15:38:48 -02:00
#if defined MHD && defined GLM
Add magnetic terms to advance_solution() and update_solution().
2011-05-01 09:44:19 -03:00
real :: decay, ch2
EVOLUTION: implement the 3rd order RK method. - implement the 3rd order Runge-Kutta time integration; - rewrite slightly the 2nd order Runge-Kutta method;
2010-12-03 15:38:48 -02:00
#endif /* MHD & GLM */
! local arrays
!
real, dimension(nqt,im,jm,km) :: u1, du
EVOLUTION: include forcing terms in the RK3 integration.
2010-12-09 17:54:59 -02:00
#ifdef FORCE
real, dimension( 3,im,jm,km) :: f
Add injected energy to total one in adiabatic case. - if we use forcing in the case of adiabatic equation of state, the total energy of each cell must be updated by the kinetic energy injected to this cell; otherwise the energy conservation will be violated;
2011-05-01 19:21:43 -03:00
#ifdef ADI
real, dimension( im,jm,km) :: ek, dek
#endif /* ADI */
EVOLUTION: include forcing terms in the RK3 integration.
2010-12-09 17:54:59 -02:00
#endif /* FORCE */
EVOLUTION: implement the 3rd order RK method. - implement the 3rd order Runge-Kutta time integration; - rewrite slightly the 2nd order Runge-Kutta method;
2010-12-03 15:38:48 -02:00
! parameters
!
real, parameter :: f4 = 1.0d0 / 4.0d0, f3 = 1.0d0 / 3.0d0
!
!-------------------------------------------------------------------------------
!
Revert "INTERPOLATION: pass the spacial increment to reconstruct()." This reverts commit d92ae5d3a62743b66e97c42be4aa5f31ef5993d1.
2010-12-06 16:20:49 -02:00
! prepare dxi, dyi, and dzi
!
dxi = adxi(pblock%meta%level)
dyi = adyi(pblock%meta%level)
dzi = adzi(pblock%meta%level)
EVOLUTION: implement the 3rd order RK method. - implement the 3rd order Runge-Kutta time integration; - rewrite slightly the 2nd order Runge-Kutta method;
2010-12-03 15:38:48 -02:00
!! 1st step of integration
!!
Revert "INTERPOLATION: pass the spacial increment to reconstruct()." This reverts commit d92ae5d3a62743b66e97c42be4aa5f31ef5993d1.
2010-12-06 16:20:49 -02:00
call update(pblock%u(:,:,:,:), du(:,:,:,:), dxi, dyi, dzi)
EVOLUTION: implement the 3rd order RK method. - implement the 3rd order Runge-Kutta time integration; - rewrite slightly the 2nd order Runge-Kutta method;
2010-12-03 15:38:48 -02:00
#ifdef SHAPE
! restrict update in a defined shape
!
call update_shapes(pblock, du(:,:,:,:))
#endif /* SHAPE */
! update the solution for the fluid variables
!
u1(1:nfl,:,:,:) = pblock%u(1:nfl,:,:,:) + dt * du(1:nfl,:,:,:)
#ifdef MHD
! update the solution for the magnetic variables
!
u1(ibx:ibz,:,:,:) = pblock%u(ibx:ibz,:,:,:) + dt * du(ibx:ibz,:,:,:)
#ifdef GLM
FORCE: update momenta due to forcinf only once. - do not put forcing terms in the RK substeps; instead, update the momenta components only once before executing the Riemann solver and update step;
2011-03-06 23:37:51 -03:00
! calculate c_h^2
!
ch2 = cmax * cmax
EVOLUTION: implement the 3rd order RK method. - implement the 3rd order Runge-Kutta time integration; - rewrite slightly the 2nd order Runge-Kutta method;
2010-12-03 15:38:48 -02:00
! update the solution for the scalar potential Psi
!
u1(iph,:,:,:) = pblock%u(iph,:,:,:) + ch2 * dt * du(iph,:,:,:)
#endif /* GLM */
#endif /* MHD */
!! 2nd step of integration
!!
Revert "INTERPOLATION: pass the spacial increment to reconstruct()." This reverts commit d92ae5d3a62743b66e97c42be4aa5f31ef5993d1.
2010-12-06 16:20:49 -02:00
call update(u1(:,:,:,:), du(:,:,:,:), dxi, dyi, dzi)
EVOLUTION: implement the 3rd order RK method. - implement the 3rd order Runge-Kutta time integration; - rewrite slightly the 2nd order Runge-Kutta method;
2010-12-03 15:38:48 -02:00
#ifdef SHAPE
! restrict update in a defined shape
!
call update_shapes(pblock, du(:,:,:,:))
#endif /* SHAPE */
! update the solution for the fluid variables
!
u1(1:nfl,:,:,:) = f4 * (3.0d0 * pblock%u(1:nfl,:,:,:) &
+ u1(1:nfl,:,:,:) + dt * du(1:nfl,:,:,:))
#ifdef MHD
! update the solution for the magnetic variables
!
u1(ibx:ibz,:,:,:) = f4 * (3.0d0 * pblock%u(ibx:ibz,:,:,:) &
+ u1(ibx:ibz,:,:,:) + dt * du(ibx:ibz,:,:,:))
#ifdef GLM
! update the solution for the scalar potential Psi
!
u1(iph,:,:,:) = f4 * (3.0d0 * pblock%u(iph,:,:,:) &
+ u1(iph,:,:,:) + ch2 * dt * du(iph,:,:,:))
#endif /* GLM */
#endif /* MHD */
!! 3rd step of integration
!!
Revert "INTERPOLATION: pass the spacial increment to reconstruct()." This reverts commit d92ae5d3a62743b66e97c42be4aa5f31ef5993d1.
2010-12-06 16:20:49 -02:00
call update(u1(:,:,:,:), du(:,:,:,:), dxi, dyi, dzi)
EVOLUTION: implement the 3rd order RK method. - implement the 3rd order Runge-Kutta time integration; - rewrite slightly the 2nd order Runge-Kutta method;
2010-12-03 15:38:48 -02:00
#ifdef SHAPE
! restrict update in a defined shape
!
call update_shapes(pblock, du(:,:,:,:))
#endif /* SHAPE */
! update the solution for the fluid variables
!
pblock%u(1:nfl,:,:,:) = f3 * (pblock%u(1:nfl,:,:,:) &
+ 2.0d0 * (u1(1:nfl,:,:,:) + dt * du(1:nfl,:,:,:)))
#ifdef MHD
! update the solution for the magnetic variables
!
pblock%u(ibx:ibz,:,:,:) = f3 * (pblock%u(ibx:ibz,:,:,:) &
+ 2.0d0 * (u1(ibx:ibz,:,:,:) + dt * du(ibx:ibz,:,:,:)))
#ifdef GLM
! update the solution for the scalar potential Psi
!
pblock%u(iph,:,:,:) = f3 * (pblock%u(iph,:,:,:) &
+ 2.0d0 * (u1(iph,:,:,:) + ch2 * dt * du(iph,:,:,:)))
! evolve analytically Psi due to the source term
!
decay = exp(- alpha_p * cmax * dt / dx_min)
pblock%u(iph,:,:,:) = decay * pblock%u(iph,:,:,:)
#endif /* GLM */
#endif /* MHD */
Add injected energy to total one in adiabatic case. - if we use forcing in the case of adiabatic equation of state, the total energy of each cell must be updated by the kinetic energy injected to this cell; otherwise the energy conservation will be violated;
2011-05-01 19:21:43 -03:00
Add forcing terms after the main update. - change order of calculating terms; first perform the main update of the conserved variables, then add forcing and source terms; in this way we do not change time step which determines the stability of scheme before the conserved variables enter Riemann solver; this is especially important in highly supersonic simulations, when in the presence of very low density, even small change of velocity can make the scheme numerically unstable;
2011-04-26 10:36:16 -03:00
#ifdef FORCE
! obtain the forcing terms in real space
!
call real_forcing(pblock%meta%level, pblock%meta%xmin, pblock%meta%ymin &
, pblock%meta%zmin, f(:,:,:,:))
Add injected energy to total one in adiabatic case. - if we use forcing in the case of adiabatic equation of state, the total energy of each cell must be updated by the kinetic energy injected to this cell; otherwise the energy conservation will be violated;
2011-05-01 19:21:43 -03:00
#ifdef ADI
! calculate kinetic energy before adding the forcing term
!
ek(:,:,:) = 0.5d0 * (pblock%u(imx,:,:,:)**2 + pblock%u(imy,:,:,:)**2 &
+ pblock%u(imz,:,:,:)**2) / pblock%u(idn,:,:,:)
#endif /* ADI */
Add forcing terms after the main update. - change order of calculating terms; first perform the main update of the conserved variables, then add forcing and source terms; in this way we do not change time step which determines the stability of scheme before the conserved variables enter Riemann solver; this is especially important in highly supersonic simulations, when in the presence of very low density, even small change of velocity can make the scheme numerically unstable;
2011-04-26 10:36:16 -03:00
! 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,:,:,:)
Add injected energy to total one in adiabatic case. - if we use forcing in the case of adiabatic equation of state, the total energy of each cell must be updated by the kinetic energy injected to this cell; otherwise the energy conservation will be violated;
2011-05-01 19:21:43 -03:00
#ifdef ADI
! calculate kinetic energy after adding the forcing term
!
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(:,:,:)
#endif /* ADI */
Add forcing terms after the main update. - change order of calculating terms; first perform the main update of the conserved variables, then add forcing and source terms; in this way we do not change time step which determines the stability of scheme before the conserved variables enter Riemann solver; this is especially important in highly supersonic simulations, when in the presence of very low density, even small change of velocity can make the scheme numerically unstable;
2011-04-26 10:36:16 -03:00
#endif /* FORCE */
EVOLUTION: implement the 3rd order RK method. - implement the 3rd order Runge-Kutta time integration; - rewrite slightly the 2nd order Runge-Kutta method;
2010-12-03 15:38:48 -02:00
!
!-------------------------------------------------------------------------------
!
end subroutine evolve_rk3
#endif /* RK3 */
Add a subroutine to calculate numerical fluxes. - the subroutine update_flux() in the SCHEME module calculates numerical fluxes from the conserved variables; the calculated fluxes are then returned;
2011-04-30 12:28:02 -03:00
#endif /* CONSERVATIVE */
Implement FLUXCT integration of the induction equation. SCHEME - implement Flux-CT scheme for the staggered magnetic field integration; BLOCK STRUCTURE - use more space efficient storage of the variables, which means storing only staggered components of magnetic field; cell-centered components are calculated only when necessary; EVOLUTION - remove loops in the field updates; operations are performed on the arrays; BOUNDARY CONDITIONS - remove loops in the bnd_copy; operations are calculated on the whole array now; INTERPOLATION - subroutine magtocen() has been rewritten to avoid problems with the array allocation; now as an argument we enter the array of all variables; subroutine uses indices for the face-centered and cell-centered magnetic field components internally; MAKE - add flag defining Flux-CT scheme; PROBLEM - use predefined array variables instead of allocated;
2010-02-28 18:35:57 -03:00
!
Updated comments in the evolution module.
2008-12-08 21:07:10 -06:00
!===============================================================================
New module 'evolution' for time integration. A new module for the time integration has been added. This module contains a set of subroutines to perform one step time integration of each leaf block using 2nd order Runge-Kutta method. More methods can be added later. Time 't', timestep 'dt' and iteration 'n' have been moved to this module as well.
2008-12-07 18:57:08 -06:00
!
end module
Reference in New Issue Copy Permalink