diff --git a/src/equations.F90 b/src/equations.F90 new file mode 100644 index 0000000..bee1b63 --- /dev/null +++ b/src/equations.F90 @@ -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 +!! +!! 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: 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 diff --git a/src/makefile b/src/makefile index 462bf46..22c54bd 100644 --- a/src/makefile +++ b/src/makefile @@ -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 \ diff --git a/src/variables.F90 b/src/variables.F90 index c12febf..d06b2e4 100644 --- a/src/variables.F90 +++ b/src/variables.F90 @@ -54,6 +54,7 @@ module variables #endif /* GLM */ #endif /* MHD */ integer(kind=4), parameter :: nvr = nqt + integer(kind=4), parameter :: nt = nqt ! !=============================================================================== !