Add new module EQUATIONS.

This module handles variable conversion for a given set of equations.
This commit is contained in:
Grzegorz Kowal 2012-07-27 16:18:02 -03:00
parent e1b79a34d0
commit 4fd6ed0f84
3 changed files with 558 additions and 6 deletions

549
src/equations.F90 Normal file
View File

@ -0,0 +1,549 @@
!!******************************************************************************
!!
!! 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-2012 Grzegorz Kowal <grzegorz@amuncode.org>
!!
!! 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 <http://www.gnu.org/licenses/>.
!!
!!******************************************************************************
!!
!! module: EQUATIONS
!!
!! This module handles supported sets of equations, provides subroutines to
!! convert between conserved and primitive variables, calculate flux and
!! characteristic speeds.
!!
!!******************************************************************************
!
module equations
! module variables are not implicit by default
!
implicit none
#ifdef ADI
! adiabatic heat ratio
!
real, save :: gamma = 5.0d0 / 3.0d0
! additional adiabatic parameters
!
real, save :: gammam1 = 2.0d0 / 3.0d0, gammam1i = 1.5d0
#endif /* ADI */
#ifdef ISO
! isothermal speed of sound and its second power
!
real, save :: csnd = 1.0d0, csnd2 = 1.0d0
#endif /* ISO */
! by default everything is public
!
public
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
contains
!
!===============================================================================
!
! subroutine INITIALIZE_EQUATIONS:
! -------------------------------
!
! Subroutine sets the default values of module parameters and obtains their
! values from the PARAMETERS module.
!
!===============================================================================
!
subroutine initialize_equations()
! include external procedures and variables
!
use parameters, only : get_parameter_real
! local variables are not implicit by default
!
implicit none
!
!-------------------------------------------------------------------------------
!
#ifdef ADI
! obtain the adiabatic specific heat ratio
!
call get_parameter_real("gamma" , gamma )
! calculate additional parameters
!
gammam1 = gamma - 1.0d0
gammam1i = 1.0d0 / gammam1
#endif /* ADI */
#ifdef ISO
! obtain the isothermal sound speed
!
call get_parameter_real("csnd" , csnd )
! calculate additional parameters
!
csnd2 = csnd * csnd
#endif /* ISO */
!-------------------------------------------------------------------------------
!
end subroutine initialize_equations
#ifdef HYDRO
!
!===============================================================================
!
! subroutine CONS2PRIM:
! --------------------
!
! Subroutine converts the conservative representation of variables to
! its corresponding primitive representation.
!
! Arguments:
!
! n - the length of input and output vectors;
! u - the input array of conservative variables;
! q - the output array of primitive variables;
!
!===============================================================================
!
subroutine cons2prim(n, u, q)
! include external procedures and variables
!
use variables , only : nt
use variables , only : idn, imx, imy, imz, ivx, ivy, ivz
#ifdef ADI
use variables , only : ipr, ien
#endif /* ADI */
! local variables are not implicit by default
!
implicit none
! input/output arguments
!
integer , intent(in) :: n
real, dimension(nt,n), intent(in) :: u
real, dimension(nt,n), intent(out) :: q
! local variables
!
integer :: i
#ifdef ADI
real :: ek, ei
#endif /* ADI */
!
!-------------------------------------------------------------------------------
!
do i = 1, n
q(idn,i) = u(idn,i)
q(ivx,i) = u(imx,i) / u(idn,i)
q(ivy,i) = u(imy,i) / u(idn,i)
q(ivz,i) = u(imz,i) / u(idn,i)
#ifdef ADI
ek = 0.5d0 * sum(u(imx:imz,i) * q(ivx:ivz,i))
ei = u(ien,i) - ek
q(ipr,i) = gammam1 * ei
#endif /* ADI */
end do
!-------------------------------------------------------------------------------
!
end subroutine cons2prim
!
!===============================================================================
!
! subroutine PRIM2CONS:
! --------------------
!
! Subroutine converts the primitive variable representation to its
! corresponding conservative representation.
!
! Arguments:
!
! n - the length of input and output vectors;
! q - the input array of primitive variables;
! u - the output array of conservative variables;
!
!===============================================================================
!
subroutine prim2cons(n, q, u)
! include external procedures and variables
!
use variables , only : nt
use variables , only : idn, imx, imy, imz, ivx, ivy, ivz
#ifdef ADI
use variables , only : ipr, ien
#endif /* ADI */
! local variables are not implicit by default
!
implicit none
! input/output arguments
!
integer , intent(in) :: n
real, dimension(nt,n), intent(in) :: q
real, dimension(nt,n), intent(out) :: u
! local variables
!
integer :: i
#ifdef ADI
real :: ek, ei
#endif /* ADI */
!
!-------------------------------------------------------------------------------
!
do i = 1, n
u(idn,i) = q(idn,i)
u(imx,i) = q(idn,i) * q(ivx,i)
u(imy,i) = q(idn,i) * q(ivy,i)
u(imz,i) = q(idn,i) * q(ivz,i)
#ifdef ADI
ek = 0.5d0 * sum(u(imx:imz,i) * q(ivx:ivz,i))
ei = gammam1i * q(ipr,i)
u(ien,i) = ei + ek
#endif /* ADI */
end do
!-------------------------------------------------------------------------------
!
end subroutine prim2cons
!
!===============================================================================
!
! subroutine FLUXSPEED:
! --------------------
!
! Subroutine calculates fluxes and characteristic speeds from primitive and
! conservative variables.
!
! Arguments:
!
! n - the length of input and output vectors;
! q - the input array of primitive variables;
! u - the input array of conservative variables;
! f - the output vector of fluxes;
! c - the output vector of characteristic speeds;
!
!===============================================================================
!
subroutine fluxspeed(n, q, u, f, c)
! include external procedures and variables
!
use variables , only : nt
use variables , only : idn, imx, imy, imz, ivx, ivy, ivz
#ifdef ADI
use variables , only : ipr, ien
#endif /* ADI */
! local variables are not implicit by default
!
implicit none
! input/output arguments
!
integer , intent(in) :: n
real, dimension(nt,n), intent(in) :: q, u
real, dimension(nt,n), intent(out) :: f
real, dimension(n) , intent(out) :: c
! local variables
!
integer :: i
!
!-------------------------------------------------------------------------------
!
do i = 1, n
! calculate the hydrodynamic fluxes
!
f(idn,i) = u(imx,i)
f(imx,i) = q(ivx,i) * u(imx,i)
f(imy,i) = q(ivx,i) * u(imy,i)
f(imz,i) = q(ivx,i) * u(imz,i)
#ifdef ISO
f(imx,i) = f(imx,i) + csnd2 * q(idn,i)
#endif /* ISO */
#ifdef ADI
f(imx,i) = f(imx,i) + q(ipr,i)
f(ien,i) = q(ivx,i) * (u(ien,i) + q(ipr,i))
#endif /* ADI */
! calculate the speed of sound
!
#ifdef ADI
c(i) = sqrt(gamma * q(ipr,i) / q(idn,i))
#endif /* ADI */
#ifdef ISO
c(i) = csnd
#endif /* ISO */
end do
!-------------------------------------------------------------------------------
!
end subroutine fluxspeed
#endif /* HYDRO */
#ifdef MHD
!
!===============================================================================
!
! subroutine CONS2PRIM:
! --------------------
!
! Subroutine converts the conservative representation of variables to
! its corresponding primitive representation.
!
! Arguments:
!
! n - the length of input and output vectors;
! u - the input array of conservative variables;
! q - the output array of primitive variables;
!
!===============================================================================
!
subroutine cons2prim(n, u, q)
! include external procedures and variables
!
use variables , only : nt
use variables , only : idn, imx, imy, imz, ivx, ivy, ivz, ibx, iby, ibz
#ifdef ADI
use variables , only : ipr, ien
#endif /* ADI */
! local variables are not implicit by default
!
implicit none
! input/output arguments
!
integer , intent(in) :: n
real, dimension(nt,n), intent(in) :: u
real, dimension(nt,n), intent(out) :: q
! local variables
!
integer :: i
#ifdef ADI
real :: ei, ek, em
#endif /* ADI */
!
!-------------------------------------------------------------------------------
!
do i = 1, n
q(idn,i) = u(idn,i)
q(ivx,i) = u(imx,i) / u(idn,i)
q(ivy,i) = u(imy,i) / u(idn,i)
q(ivz,i) = u(imz,i) / u(idn,i)
q(ibx,i) = u(ibx,i)
q(iby,i) = u(iby,i)
q(ibz,i) = u(ibz,i)
#ifdef ADI
ek = 0.5d0 * sum(u(imx:imz,i) * q(ivx:ivz,i))
em = 0.5d0 * sum(q(ibx:ibz,i) * q(ibx:ibz,i))
ei = u(ien,i) - (ek + em)
q(ipr,i) = gammam1 * ei
#endif /* ADI */
end do
!-------------------------------------------------------------------------------
!
end subroutine cons2prim
!
!===============================================================================
!
! subroutine PRIM2CONS:
! --------------------
!
! Subroutine converts the primitive variable representation to its
! corresponding conservative representation.
!
! Arguments:
!
! n - the length of input and output vectors;
! q - the input array of primitive variables;
! u - the output array of conservative variables;
!
!===============================================================================
!
subroutine prim2cons(n, q, u)
! include external procedures and variables
!
use variables , only : nt
use variables , only : idn, imx, imy, imz, ivx, ivy, ivz, ibx, iby, ibz
#ifdef ADI
use variables , only : ipr, ien
#endif /* ADI */
! local variables are not implicit by default
!
implicit none
! input/output arguments
!
integer , intent(in) :: n
real, dimension(nt,n), intent(in) :: q
real, dimension(nt,n), intent(out) :: u
! local variables
!
integer :: i
#ifdef ADI
real :: ei, ek, em
#endif /* ADI */
!
!-------------------------------------------------------------------------------
!
do i = 1, n
u(idn,i) = q(idn,i)
u(imx,i) = q(idn,i) * q(ivx,i)
u(imy,i) = q(idn,i) * q(ivy,i)
u(imz,i) = q(idn,i) * q(ivz,i)
u(ibx,i) = q(ibx,i)
u(iby,i) = q(iby,i)
u(ibz,i) = q(ibz,i)
#ifdef ADI
ei = gammam1i * q(ipr,i)
ek = 0.5d0 * sum(u(imx:imz,i) * q(ivx:ivz,i))
em = 0.5d0 * sum(q(ibx:ibz,i) * q(ibx:ibz,i))
u(ien,i) = ei + ek + em
#endif /* ADI */
end do
!-------------------------------------------------------------------------------
!
end subroutine prim2cons
!
!===============================================================================
!
! subroutine FLUXSPEED:
! --------------------
!
! Subroutine calculates fluxes and characteristic speeds from primitive and
! conservative variables.
!
! Arguments:
!
! n - the length of input and output vectors;
! q - the input array of primitive variables;
! u - the input array of conservative variables;
! f - the output vector of fluxes;
! c - the output vector of characteristic speeds;
!
!===============================================================================
!
subroutine fluxspeed(n, q, u, f, c)
! include external procedures and variables
!
use variables , only : nt
use variables , only : idn, imx, imy, imz, ivx, ivy, ivz
use variables , only : ibx, iby, ibz
#ifdef ADI
use variables , only : ipr, ien
#endif /* ADI */
! local variables are not implicit by default
!
implicit none
! input/output arguments
!
integer , intent(in) :: n
real, dimension(nt,n), intent(in) :: q, u
real, dimension(nt,n), intent(out) :: f
real, dimension(n) , intent(out) :: c
! local variables
!
integer :: i
real :: bb, vb, pm, pt
real :: cs2, ca2, cx2, cf2, cl2
!
!-------------------------------------------------------------------------------
!
do i = 1, n
! prepare pressures and scalar product
!
bb = sum(q(ibx:ibz,i) * q(ibx:ibz,i))
pm = 0.5d0 * bb
#ifdef ADI
vb = sum(q(ivx:ivz,i) * q(ibx:ibz,i))
pt = q(ipr,i)
#endif /* ADI */
#ifdef ISO
pt = csnd2 * q(idn,i)
#endif /* ISO */
pt = pt + pm
! calculate the magnetohydrodynamic fluxes
!
f(idn,i) = u(imx,i)
f(imx,i) = q(ivx,i) * u(imx,i) - q(ibx,i) * q(ibx,i)
f(imy,i) = q(ivx,i) * u(imy,i) - q(ibx,i) * q(iby,i)
f(imz,i) = q(ivx,i) * u(imz,i) - q(ibx,i) * q(ibz,i)
f(imx,i) = f(imx,i) + pt
#ifdef ADI
f(ien,i) = q(ivx,i) * (u(ien,i) + pt) - q(ibx,i) * vb
#endif /* ADI */
f(ibx,i) = 0.0d0
f(iby,i) = q(ivx,i) * q(iby,i) - q(ibx,i) * q(ivy,i)
f(ibz,i) = q(ivx,i) * q(ibz,i) - q(ibx,i) * q(ivz,i)
! calculate the fast magnetosonic speed
!
#ifdef ADI
cs2 = gamma * q(ipr,i) / q(idn,i)
#endif /* ADI */
#ifdef ISO
cs2 = csnd2
#endif /* ISO */
ca2 = sum(q(ibx:ibz,i) * q(ibx:ibz,i)) / q(idn,i)
cx2 = q(ibx,i) * q(ibx,i) / q(idn,i)
cf2 = cs2 + ca2
cl2 = max(0.0d0, cf2**2 - 4.0d0 * cs2 * cx2)
c(i) = sqrt(0.5d0 * (cf2 + sqrt(cl2)))
end do
!-------------------------------------------------------------------------------
!
end subroutine fluxspeed
#endif /* MHD */
!===============================================================================
!
end module equations

View File

@ -161,13 +161,14 @@ name = amun
default: $(name).x
sources = blocks.F90 boundaries.F90 config.F90 constants.F90 coordinates.F90 \
driver.F90 error.F90 evolution.F90 forcing.F90 integrals.F90 \
interpolation.F90 io.F90 mesh.F90 mpitools.F90 parameters.F90 \
problem.F90 random.F90 scheme.F90 timers.F90 variables.F90
driver.F90 equations.F90 error.F90 evolution.F90 forcing.F90 \
integrals.F90 interpolation.F90 io.F90 mesh.F90 mpitools.F90 \
parameters.F90 problem.F90 random.F90 scheme.F90 timers.F90 \
variables.F90
objects = blocks.o boundaries.o config.o constants.o coordinates.o driver.o \
error.o evolution.o forcing.o integrals.o interpolation.o io.o \
mesh.o mpitools.o parameters.o problem.o random.o scheme.o timers.o \
variables.o
equations.o error.o evolution.o forcing.o integrals.o \
interpolation.o io.o mesh.o mpitools.o parameters.o problem.o \
random.o scheme.o timers.o variables.o
files = $(sources) makefile make.default config.in license.txt hosts
$(name).x: $(objects)
@ -210,6 +211,7 @@ coordinates.o : coordinates.F90 parameters.o
driver.o : driver.F90 blocks.o config.o coordinates.o evolution.o \
forcing.o integrals.o io.o mesh.o mpitools.o parameters.o \
random.o
equations.o : equations.F90 parameters.o variables.o
error.o : error.F90
evolution.o : evolution.F90 blocks.o boundaries.o config.o coordinates.o \
forcing.o interpolation.o mesh.o mpitools.o problem.o \

View File

@ -54,6 +54,7 @@ module variables
#endif /* GLM */
#endif /* MHD */
integer(kind=4), parameter :: nvr = nqt
integer(kind=4), parameter :: nt = nqt
!
!===============================================================================
!