!!****************************************************************************** !! !! This file is part of the AMUN source code, a program to perform !! Newtonian or relativistic magnetohydrodynamical simulations on uniform or !! adaptive mesh. !! !! Copyright (C) 2008-2014 Grzegorz Kowal !! !! 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 3 of the License, or !! (at your option) any later version. !! !! This program is distributed in the hope that it will be useful, !! 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 !! along with this program. If not, see . !! !!****************************************************************************** !! !! module: EVOLUTION !! !! This module provides an interface for temporal integration with !! the stability handling. New integration methods can be added by !! implementing more evolve_* subroutines. !! !!****************************************************************************** ! module evolution ! module variables are not implicit by default ! implicit none ! pointer to the temporal integration subroutine ! procedure(evolve_euler), pointer, save :: evolve => null() ! evolution parameters ! real , save :: cfl = 0.5d+00 ! coefficient controlling the decay of scalar potential ѱ ! real , save :: alpha = 2.0d+00 real , save :: decay = 1.0d+00 ! time variables ! integer, save :: step = 0 real , save :: time = 0.0d+00 real , save :: dt = 1.0d+00 real , save :: dtn = 1.0d+00 ! by default everything is private ! private ! declare public subroutines ! public :: initialize_evolution, finalize_evolution public :: advance, new_time_step ! declare public variables ! public :: cfl, step, time, dt, dtn !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! contains ! !=============================================================================== !! !!*** PUBLIC SUBROUTINES ***************************************************** !! !=============================================================================== ! !=============================================================================== ! ! subroutine INITIALIZE_EVOLUTION: ! ------------------------------- ! ! Subroutine initializes module EVOLUTION by setting its parameters. ! ! Arguments: ! ! verbose - a logical flag turning the information printing; ! iret - an integer flag for error return value; ! !=============================================================================== ! subroutine initialize_evolution(verbose, iret) ! include external procedures ! use parameters, only : get_parameter_string, get_parameter_real ! local variables are not implicit by default ! implicit none ! subroutine arguments ! logical, intent(in) :: verbose integer, intent(inout) :: iret ! local variables ! character(len=255) :: integration = "rk2" character(len=255) :: name_int = "" ! !------------------------------------------------------------------------------- ! ! get the integration method and the value of the CFL coefficient ! call get_parameter_string("time_advance", integration) call get_parameter_real ("cfl" , cfl ) call get_parameter_real ("alpha" , alpha ) ! select the integration method, check the correctness of the integration ! parameters and adjust the CFL coefficient if necessary ! select case(trim(integration)) case ("euler", "EULER") name_int = "1st order Euler method" evolve => evolve_euler case ("rk2", "RK2") name_int = "2nd order Runge-Kutta method" evolve => evolve_rk2 case default if (verbose) then write (*,"(1x,a)") "The selected time advance method is not " // & "implemented: " // trim(integration) name_int = "2nd order Runge-Kutta method" evolve => evolve_rk2 end if end select ! calculate the decay factor for magnetic field divergence scalar source term ! decay = exp(- alpha * cfl) ! print information about the Riemann solver ! if (verbose) then write (*,"(4x,a,1x,a)" ) "time advance =", trim(name_int) end if !------------------------------------------------------------------------------- ! end subroutine initialize_evolution ! !=============================================================================== ! ! subroutine FINALIZE_EVOLUTION: ! ----------------------------- ! ! Subroutine releases memory used by the module. ! ! Arguments: ! ! iret - an integer flag for error return value; ! !=============================================================================== ! subroutine finalize_evolution(iret) ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer, intent(inout) :: iret ! !------------------------------------------------------------------------------- ! nullify(evolve) !------------------------------------------------------------------------------- ! end subroutine finalize_evolution ! !=============================================================================== ! ! subroutine ADVANCE: ! ------------------ ! ! Subroutine advances the solution by one time step using the selected time ! integration method. ! ! !=============================================================================== ! subroutine advance() ! include external procedures ! use boundaries , only : boundary_variables use mesh , only : update_mesh ! include external variables ! use coordinates , only : toplev ! local variables are not implicit by default ! implicit none ! !------------------------------------------------------------------------------- ! ! find new time step ! call new_time_step() ! advance the solution using the selected method ! call evolve() ! chec if we need to perform the refinement step ! if (toplev > 1) then ! check refinement and refine ! call update_mesh() ! update boundaries ! call boundary_variables() end if ! toplev > 1 ! update primitive variables ! call update_variables() !------------------------------------------------------------------------------- ! end subroutine advance ! !=============================================================================== !! !!*** PRIVATE SUBROUTINES **************************************************** !! !=============================================================================== ! !=============================================================================== ! ! subroutine EVOLVE_EULER: ! ----------------------- ! ! Subroutine advances the solution by one time step using the 1st order ! Euler integration method. ! !=============================================================================== ! subroutine evolve_euler() ! include external procedures ! use boundaries , only : boundary_variables use schemes , only : update_increment ! include external variables ! use blocks , only : block_data, list_data use coordinates , only : adx, ady, adz use coordinates , only : im, jm, km use equations , only : nv, ibp ! local variables are not implicit by default ! implicit none ! local pointers ! type(block_data), pointer :: pblock ! local arrays ! real, dimension(3) :: dh real, dimension(nv,im,jm,km) :: du ! !------------------------------------------------------------------------------- ! ! update fluxes for the first step of the RK2 integration ! call update_fluxes() ! update the solution using numerical fluxes stored in the data blocks ! pblock => list_data do while (associated(pblock)) ! obtain dx, dy, and dz for the current block ! dh(1) = dt / adx(pblock%meta%level) dh(2) = dt / ady(pblock%meta%level) dh(3) = dt / adz(pblock%meta%level) ! calculate variable increment for the current block ! call update_increment(dh(:), pblock%f(:,:,:,:,:), du(:,:,:,:)) ! update the solution for the fluid variables ! pblock%u0(1:nv,:,:,:) = pblock%u0(1:nv,:,:,:) + du(1:nv,:,:,:) ! update the conservative variable pointer ! pblock%u => pblock%u0 ! update ψ by its source term ! if (ibp > 0) pblock%u(ibp,:,:,:) = decay * pblock%u(ibp,:,:,:) ! assign pointer to the next block ! pblock => pblock%next end do ! update boundaries ! call boundary_variables() ! update primitive variables ! call update_variables() !------------------------------------------------------------------------------- ! end subroutine evolve_euler ! !=============================================================================== ! ! subroutine EVOLVE_RK2: ! --------------------- ! ! Subroutine advances the solution by one time step using the 2nd order ! Runge-Kutta time integration method. ! ! !=============================================================================== ! subroutine evolve_rk2() ! include external procedures ! use boundaries , only : boundary_variables use schemes , only : update_increment ! include external variables ! use blocks , only : block_data, list_data use coordinates , only : adx, ady, adz use coordinates , only : im, jm, km use equations , only : nv, ibp ! local variables are not implicit by default ! implicit none ! local pointers ! type(block_data), pointer :: pblock ! local arrays ! real, dimension(3) :: dh real, dimension(nv,im,jm,km) :: du ! !------------------------------------------------------------------------------- ! ! update fluxes for the first step of the RK2 integration ! call update_fluxes() ! update the solution using numerical fluxes stored in the data blocks ! pblock => list_data do while (associated(pblock)) ! obtain dx, dy, and dz for the current block ! dh(1) = dt / adx(pblock%meta%level) dh(2) = dt / ady(pblock%meta%level) dh(3) = dt / adz(pblock%meta%level) ! calculate variable increment for the current block ! call update_increment(dh(:), pblock%f(:,:,:,:,:), du(:,:,:,:)) ! update the solution for the fluid variables ! pblock%u1(1:nv,:,:,:) = pblock%u0(1:nv,:,:,:) + du(1:nv,:,:,:) ! update the conservative variable pointer ! pblock%u => pblock%u1 ! assign pointer to the next block ! pblock => pblock%next end do ! update boundaries ! call boundary_variables() ! update primitive variables ! call update_variables() ! update fluxes for the second step of the RK2 integration ! call update_fluxes() ! update the solution using numerical fluxes stored in the data blocks ! pblock => list_data do while (associated(pblock)) ! obtain dx, dy, and dz for the current block ! dh(1) = dt / adx(pblock%meta%level) dh(2) = dt / ady(pblock%meta%level) dh(3) = dt / adz(pblock%meta%level) ! calculate variable increment for the current block ! call update_increment(dh(:), pblock%f(:,:,:,:,:), du(:,:,:,:)) ! update the solution for the fluid variables ! pblock%u0(1:nv,:,:,:) = 0.5d0 * (pblock%u0(1:nv,:,:,:) & + pblock%u1(1:nv,:,:,:) + du(1:nv,:,:,:)) ! update the conservative variable pointer ! pblock%u => pblock%u0 ! update ψ by its source term ! if (ibp > 0) pblock%u(ibp,:,:,:) = decay * pblock%u(ibp,:,:,:) ! assign pointer to the next block ! pblock => pblock%next end do ! update boundaries ! call boundary_variables() ! update primitive variables ! call update_variables() !------------------------------------------------------------------------------- ! end subroutine evolve_rk2 ! !=============================================================================== ! ! subroutine UPDATE_FLUXES: ! ------------------------ ! ! Subroutine iterates over all data blocks and calculates the numerical ! fluxes for each block. After the fluxes are updated, they are corrected ! for blocks which have neighbours at higher refinement level. ! ! !=============================================================================== ! subroutine update_fluxes() ! include external procedures ! use boundaries , only : boundary_fluxes use schemes , only : update_flux ! include external variables ! use blocks , only : block_data, list_data use coordinates , only : adx, ady, adz ! local variables are not implicit by default ! implicit none ! local pointers ! type(block_data), pointer :: pblock ! local vectors ! real, dimension(3) :: dx ! local variables ! integer :: n ! !------------------------------------------------------------------------------- ! ! iterate over all data blocks ! pblock => list_data do while (associated(pblock)) ! obtain dx, dy, and dz for the current block ! dx(1) = adx(pblock%meta%level) dx(2) = ady(pblock%meta%level) dx(3) = adz(pblock%meta%level) ! update the flux for the current block ! do n = 1, NDIMS call update_flux(n, dx(n), pblock%q(:,:,:,:), pblock%f(n,:,:,:,:)) end do ! assign pointer to the next block ! pblock => pblock%next end do ! correct the numerical fluxes of the blocks which have neighbours at higher ! level ! call boundary_fluxes() !------------------------------------------------------------------------------- ! end subroutine update_fluxes ! !=============================================================================== ! ! subroutine UPDATE_VARIABLES: ! --------------------------- ! ! Subroutine iterates over all data blocks and converts the conservative ! variables to their primitive representation. ! ! !=============================================================================== ! subroutine update_variables() ! include external procedures ! use equations , only : update_primitive_variables ! include external variables ! use blocks , only : block_data, list_data ! local variables are not implicit by default ! implicit none ! local pointers ! type(block_data), pointer :: pblock ! !------------------------------------------------------------------------------- ! ! iterate over all data blocks ! pblock => list_data do while (associated(pblock)) ! convert conserved variables to primitive ones for the current block ! call update_primitive_variables(pblock%u, pblock%q) ! assign pointer to the next block ! pblock => pblock%next end do !------------------------------------------------------------------------------- ! end subroutine update_variables ! !=============================================================================== ! ! subroutine NEW_TIME_STEP: ! ------------------------ ! ! Subroutine estimates the new time step from the maximum speed in the system ! and source term constraints. ! ! !=============================================================================== ! subroutine new_time_step() ! include external procedures ! use equations , only : maxspeed, cmax, cmax2 #ifdef MPI use mpitools , only : reduce_maximum_real, reduce_maximum_integer #endif /* MPI */ ! include external variables ! use blocks , only : block_data, list_data use coordinates , only : adx, ady, adz use coordinates , only : toplev ! local variables are not implicit by default ! implicit none ! local pointers ! type(block_data), pointer :: pblock ! local variables ! integer :: iret integer(kind=4) :: lev real :: cm, dx_min ! local parameters ! real, parameter :: eps = tiny(cmax) ! !------------------------------------------------------------------------------- ! ! reset the maximum speed, and the highest level ! cmax = eps lev = 1 ! iterate over all data blocks in order to find the maximum speed among them ! and the highest level which is required to obtain the spatial step ! pblock => list_data do while (associated(pblock)) ! find the maximum level occupied by blocks (can be smaller than toplev) ! lev = max(lev, pblock%meta%level) ! obtain the maximum speed for the current block ! cm = maxspeed(pblock%q(:,:,:,:)) ! compare global and local maximum speeds ! cmax = max(cmax, cm) ! assiociate the pointer with the next block ! pblock => pblock%next end do #ifdef MPI ! find maximum speed in the system from all processors ! call reduce_maximum_real (cmax, iret) call reduce_maximum_integer(lev , iret) #endif /* MPI */ ! calculate squared cmax ! cmax2 = cmax * cmax ! find the smallest spatial step ! #if NDIMS == 2 dx_min = min(adx(lev), ady(lev)) #endif /* NDIMS == 2 */ #if NDIMS == 3 dx_min = min(adx(lev), ady(lev), adz(lev)) #endif /* NDIMS == 3 */ ! calcilate the new time step ! dtn = dx_min / max(cmax, eps) ! calculate the new time step ! dt = cfl * dtn !------------------------------------------------------------------------------- ! end subroutine new_time_step !=============================================================================== ! end module evolution