!!****************************************************************************** !! !! 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-2024 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: PROBLEMS !! !! This module handles the initialization of various test and research !! problems. !! !! !!****************************************************************************** ! module problems implicit none character(len=64), save :: problem_name = "none" private public :: initialize_problems, finalize_problems public :: problem_name !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! contains ! !=============================================================================== !! !!*** PUBLIC SUBROUTINES ***************************************************** !! !=============================================================================== ! !=============================================================================== ! ! subroutine INITIALIZE_PROBLEMS: ! ------------------------------ ! ! Subroutine prepares module PROBLEMS. ! ! Arguments: ! ! problem - the problem name ! rcount - the run count for restarted jobs ! verbose - a logical flag turning the information printing; ! status - an integer flag for error return value; ! !=============================================================================== ! subroutine initialize_problems(problem, rcount, verbose, status) use mesh , only : setup_problem use parameters , only : get_parameter implicit none character(len=64), intent(in) :: problem integer , intent(in) :: rcount logical , intent(in) :: verbose integer , intent(out) :: status !------------------------------------------------------------------------------- ! status = 0 problem_name = problem select case(trim(problem)) case("riemann") setup_problem => setup_problem_riemann case("blast") setup_problem => setup_problem_blast case("st", "sedov-taylor", "ST", "Sedov-Taylor") setup_problem => setup_problem_sedov_taylor case("implosion") setup_problem => setup_problem_implosion case("kh", "kelvinhelmholtz", "kelvin-helmholtz", "Kelvin-Helmholtz") setup_problem => setup_problem_kelvin_helmholtz case("rt", "rayleightaylor", "rayleigh-taylor") setup_problem => setup_problem_rayleigh_taylor case("ot", "orszag-tang", "Orszag-Tang") setup_problem => setup_problem_orszag_tang case("current-sheet", "current_sheet") setup_problem => setup_problem_current_sheet case("tearing-mode", "tearing_mode", "tearing") setup_problem => setup_problem_tearing case("turbulence") setup_problem => setup_problem_turbulence case("stellar_wind", "stellar-wind", "wind") setup_problem => setup_problem_wind case default end select !------------------------------------------------------------------------------- ! end subroutine initialize_problems ! !=============================================================================== ! ! subroutine FINALIZE_PROBLEMS: ! ---------------------------- ! ! Subroutine releases memory used by the module. ! ! Arguments: ! ! status - an integer flag for error return value; ! !=============================================================================== ! subroutine finalize_problems(status) implicit none integer, intent(out) :: status !------------------------------------------------------------------------------- ! status = 0 !------------------------------------------------------------------------------- ! end subroutine finalize_problems ! !=============================================================================== !! !!*** PRIVATE SUBROUTINES **************************************************** !! !=============================================================================== ! !=============================================================================== ! ! subroutine SETUP_PROBLEM_RIEMANN: ! -------------------------------- ! ! Subroutine sets the initial conditions for the general Riemann problem. ! ! Arguments: ! ! pdata - pointer to the datablock structure of the currently initialized ! block; ! ! !=============================================================================== ! subroutine setup_problem_riemann(pdata) use blocks , only : block_data use coordinates, only : nn => bcells use coordinates, only : ax, adx use equations , only : prim2cons use equations , only : nv use equations , only : qpbnd implicit none type(block_data), pointer, intent(inout) :: pdata integer :: p, i, j, k real(kind=8) :: xl, xr real(kind=8) :: dx, dxh real(kind=8), dimension(nv,nn) :: q, u real(kind=8), dimension(nn) :: x ! !------------------------------------------------------------------------------- ! ! prepare block coordinates ! x(:) = pdata%meta%xmin + ax(pdata%meta%level,:) ! calculate mesh intervals and areas ! dx = adx(pdata%meta%level) dxh = 0.5d+00 * dx ! set the left and right states of the primitive variables ! do i = 1, nn xl = x(i) - dxh xr = x(i) + dxh if (xr <= 0.0d+00) then do p = 1, nv q(p,i) = qpbnd(p,1,1) end do else if (xl >= 0.0d+00) then do p = 1, nv q(p,i) = qpbnd(p,1,2) end do else do p = 1, nv q(p,i) = (xr * qpbnd(p,1,2) - xl * qpbnd(p,1,1)) / dx end do end if end do ! i = 1, im ! convert the primitive variables to conservative ones ! call prim2cons(q(:,:), u(:,:)) ! iterate over all positions in the YZ plane ! #if NDIMS == 2 k = 1 #endif /* NDIMS == 2 */ #if NDIMS == 3 do k = 1, nn #endif /* NDIMS == 3 */ do j = 1, nn ! copy the conserved variables to the current block ! pdata%u(:,:,j,k) = u(:,:) ! copy the primitive variables to the current block ! pdata%q(:,:,j,k) = q(:,:) end do ! j = 1, nn #if NDIMS == 3 end do ! k = 1, nn #endif /* NDIMS == 3 */ !------------------------------------------------------------------------------- ! end subroutine setup_problem_riemann ! !=============================================================================== ! ! subroutine SETUP_PROBLEM_BLAST: ! ------------------------------ ! ! Subroutine sets the initial conditions for the spherical blast problem. ! ! Arguments: ! ! pdata - pointer to the datablock structure of the currently initialized ! block; ! ! !=============================================================================== ! subroutine setup_problem_blast(pdata) ! include external procedures and variables ! use blocks , only : block_data use constants , only : d2r use coordinates, only : nn => bcells use coordinates, only : ax, ay, adx, ady, advol #if NDIMS == 3 use coordinates, only : az, adz #endif /* NDIMS == 3 */ use equations , only : prim2cons use equations , only : adiabatic_index use equations , only : nv use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp use parameters , only : get_parameter ! local variables are not implicit by default ! implicit none ! input arguments ! type(block_data), pointer, intent(inout) :: pdata ! default parameter values ! real(kind=8), save :: dens = 1.00d+00 real(kind=8), save :: ratio = 1.00d+02 real(kind=8), save :: radius = 1.00d-01 real(kind=8), save :: csnd = 4.0824829046386301635d-01 real(kind=8), save :: buni = 1.00d+00 real(kind=8), save :: angle = 4.50d+01 #if NDIMS == 3 integer , save :: nsubgrid = 10 #endif /* NDIMS == 3 */ ! local saved parameters ! logical , save :: first = .true. real(kind=8), save :: dn_amb, dn_ovr real(kind=8), save :: pr_amb, pr_ovr real(kind=8), save :: bx, by real(kind=8), save :: r2 ! local variables ! integer :: i, j, k #if NDIMS == 3 integer :: ic, jc, kc #endif /* NDIMS == 3 */ real(kind=8) :: xl, yl, xu, yu, rl, ru real(kind=8) :: dx, dy, dxh, dyh, dvol real(kind=8) :: sn #if NDIMS == 3 real(kind=8) :: zl, zu, dz, dzh real(kind=8) :: xb, yb, zb real(kind=8) :: xt, yt, zt real(kind=8) :: fc_inc #else /* NDIMS == 3 */ real(kind=8) :: rlu, rul real(kind=8) :: xb, yb real(kind=8) :: xt, yt real(kind=8) :: ph #endif /* NDIMS == 3 */ real(kind=8) :: fc_amb, fc_ovr ! local arrays ! real(kind=8), dimension(nv,nn) :: q, u real(kind=8), dimension(nn) :: x, y #if NDIMS == 3 real(kind=8), dimension(nn) :: z ! allocatable arrays ! real(kind=8), dimension(:), allocatable :: xm, ym, zm real(kind=8), dimension(:), allocatable :: xp, yp, zp #endif /* NDIMS == 3 */ ! !------------------------------------------------------------------------------- ! ! prepare problem constants during the first subroutine call ! if (first) then ! get problem parameters ! call get_parameter("dens" , dens ) call get_parameter("ratio" , ratio ) call get_parameter("radius" , radius) call get_parameter("sound_speed", csnd ) call get_parameter("buni" , buni ) call get_parameter("angle" , angle ) #if NDIMS == 3 ! get the fine grid resolution ! call get_parameter("nsubgrid", nsubgrid) ! correct subgrid resolution if necessary ! nsubgrid = max(1, nsubgrid) #endif /* NDIMS == 3 */ ! calculate the overdense and ambient region densities ! dn_amb = dens if (ipr > 0) then dn_ovr = dn_amb ! calculate parallel and perpendicular pressures from sound speeds ! pr_amb = dens * csnd * csnd / adiabatic_index pr_ovr = pr_amb * ratio else dn_ovr = dn_amb * ratio end if ! calculate initial uniform field components ! if (ibx > 0) then sn = sin(d2r * angle) bx = buni * sqrt(1.0d+00 - sn * sn) by = buni * sn end if ! calculate the square of radius ! r2 = radius * radius ! reset the first execution flag ! first = .false. end if ! first call ! prepare block coordinates ! x(:) = pdata%meta%xmin + ax(pdata%meta%level,:) y(:) = pdata%meta%ymin + ay(pdata%meta%level,:) #if NDIMS == 3 z(:) = pdata%meta%zmin + az(pdata%meta%level,:) #endif /* NDIMS == 3 */ ! calculate mesh intervals and areas ! dx = adx(pdata%meta%level) dy = ady(pdata%meta%level) #if NDIMS == 3 dz = adz(pdata%meta%level) #endif /* NDIMS == 3 */ dxh = 0.5d+00 * dx dyh = 0.5d+00 * dy #if NDIMS == 3 dzh = 0.5d+00 * dz #endif /* NDIMS == 3 */ dvol = advol(pdata%meta%level) #if NDIMS == 3 ! allocate subgrid coordinates ! allocate(xm(nsubgrid), ym(nsubgrid), zm(nsubgrid)) allocate(xp(nsubgrid), yp(nsubgrid), zp(nsubgrid)) ! and generate them ! xm(:) = (1.0d+00 * (/(i, i = 0, nsubgrid - 1)/)) / nsubgrid ym(:) = xm(:) zm(:) = xm(:) xm(:) = xm(:) * dx ym(:) = ym(:) * dy zm(:) = zm(:) * dz xp(:) = (1.0d+00 * (/(i, i = 1, nsubgrid )/)) / nsubgrid yp(:) = xp(:) zp(:) = xp(:) xp(:) = xp(:) * dx yp(:) = yp(:) * dy zp(:) = zp(:) * dz ! calculate the factor increment for the given subgrid ! fc_inc = dvol / nsubgrid**3 #endif /* NDIMS == 3 */ ! set the ambient density and pressure ! q(idn,:) = dn_amb if (ipr > 0) q(ipr,:) = pr_amb ! reset velocity components ! q(ivx,:) = 0.0d+00 q(ivy,:) = 0.0d+00 q(ivz,:) = 0.0d+00 ! if magnetic field is present, set it to be uniform with the desired strength ! and orientation ! if (ibx > 0) then q(ibx,:) = bx q(iby,:) = by q(ibz,:) = 0.0d+00 q(ibp,:) = 0.0d+00 end if ! iterate over all positions in the YZ plane ! #if NDIMS == 2 k = 1 #endif /* NDIMS == 2 */ #if NDIMS == 3 do k = 1, nn ! calculate the corner Z coordinates ! zl = abs(z(k)) - dzh zu = abs(z(k)) + dzh #endif /* NDIMS == 3 */ do j = 1, nn ! calculate the corner Y coordinates ! yl = abs(y(j)) - dyh yu = abs(y(j)) + dyh ! sweep along the X coordinate ! do i = 1, nn ! calculate the corner X coordinates ! xl = abs(x(i)) - dxh xu = abs(x(i)) + dxh ! calculate the minimum and maximum corner distances from the origin ! #if NDIMS == 3 rl = xl * xl + yl * yl + zl * zl ru = xu * xu + yu * yu + zu * zu #else /* NDIMS == 3 */ rl = xl * xl + yl * yl ru = xu * xu + yu * yu #endif /* NDIMS == 3 */ ! set the initial density and pressure in cells laying completely within ! the blast radius ! if (ru <= r2) then ! set the overpressure region density ! q(idn,i) = dn_ovr ! set the overpressure region pressure ! if (ipr > 0) q(ipr,i) = pr_ovr ! set the initial pressure in the cell completely outside the radius ! else if (rl >= r2) then ! set the ambient region density ! q(idn,i) = dn_amb ! set the ambient medium pressure ! if (ipr > 0) q(ipr,i) = pr_amb ! integrate density or pressure in cells which are crossed by the circule with ! the given radius ! else #if NDIMS == 3 ! interpolate the factor using subgrid ! fc_ovr = 0.0d+00 do kc = 1, nsubgrid zb = (zl + zm(kc))**2 zt = (zl + zp(kc))**2 do jc = 1, nsubgrid yb = (yl + ym(jc))**2 yt = (yl + yp(jc))**2 do ic = 1, nsubgrid xb = (xl + xm(ic))**2 xt = (xl + xp(ic))**2 ! update the integration factor depending on the subcell position ! if ((xt + yt + zt) <= r2) then fc_ovr = fc_ovr + fc_inc else if ((xb + yb + zb) < r2) then fc_ovr = fc_ovr + 0.5d+00 * fc_inc end if end do ! ic = 1, nsubgrid end do ! jc = 1, nsubgrid end do ! kc = 1, nsubgrid #else /* NDIMS == 3 */ ! calculate the distance of remaining corners ! rlu = xl * xl + yu * yu rul = xu * xu + yl * yl ! separate in the cases of which corners lay inside, and which outside ! the radius ! if (min(rlu, rul) >= r2) then ! only one cell corner inside the radius ! ! calculate middle coordinates of the radius-edge crossing point ! xb = sqrt(r2 - yl**2) - xl yb = sqrt(r2 - xl**2) - yl ! calculate the sin(½φ), φ, and sin(φ) ! sn = 0.5d+00 * sqrt(xb**2 + yb**2) / radius ph = 2.0d+00 * asin(sn) sn = sin(ph) ! calculate the area of cell intersection with the radius ! fc_ovr = 0.5d+00 * (xb * yb + (ph - sn) * r2) else if (rlu >= r2) then ! two lower corners inside the radius ! ! calculate middle coordinates of the radius-edge crossing point ! yb = sqrt(r2 - xl**2) - yl yt = sqrt(r2 - xu**2) - yl ! calculate the sin(½φ), φ, and sin(φ) ! sn = 0.5d+00 * sqrt(dx**2 + (yt - yb)**2) / radius ph = 2.0d+00 * asin(sn) sn = sin(ph) ! calculate the area of cell intersection with the radius ! fc_ovr = 0.5d+00 * ((yt + yb) * dx + (ph - sn) * r2) else if (rul >= r2) then ! two left corners inside the radius ! ! calculate middle coordinates of the radius-edge crossing point ! xb = sqrt(r2 - yl**2) - xl xt = sqrt(r2 - yu**2) - xl ! calculate the sin(½φ), φ, and sin(φ) ! sn = 0.5d+00 * sqrt((xt - xb)**2 + dy**2) / radius ph = 2.0d+00 * asin(sn) sn = sin(ph) ! calculate the area of cell intersection with the radius ! fc_ovr = 0.5d+00 * ((xt + xb) * dy + (ph - sn) * r2) else ! three corners inside the radius ! ! calculate middle coordinates of the radius-edge crossing point ! xt = xu - sqrt(r2 - yu**2) yt = yu - sqrt(r2 - xu**2) ! calculate the sin(½φ), φ, and sin(φ) ! sn = 0.5d+00 * sqrt(xt**2 + yt**2) / radius ph = 2.0d+00 * asin(sn) sn = sin(ph) ! calculate the area of cell intersection with the radius ! fc_ovr = dvol - 0.5d+00 * (xt * yt - (ph - sn) * r2) end if #endif /* NDIMS == 3 */ ! normalize coefficients ! fc_ovr = fc_ovr / dvol fc_amb = 1.0d+00 - fc_ovr ! integrate the density over the edge cells ! q(idn,i) = fc_ovr * dn_ovr + fc_amb * dn_amb ! integrate the pressure over the edge cells ! if (ipr > 0) q(ipr,i) = fc_ovr * pr_ovr + fc_amb * pr_amb end if end do ! i = 1, nn ! convert the primitive variables to conservative ones ! call prim2cons(q(:,:), u(:,:)) ! copy the conserved variables to the current block ! pdata%u(:,:,j,k) = u(:,:) ! copy the primitive variables to the current block ! pdata%q(:,:,j,k) = q(:,:) end do ! j = 1, nn #if NDIMS == 3 end do ! k = 1, nn ! deallocate subgrid coordinates ! deallocate(xm, ym, zm) deallocate(xp, yp, zp) #endif /* NDIMS == 3 */ !------------------------------------------------------------------------------- ! end subroutine setup_problem_blast ! !=============================================================================== ! ! subroutine SETUP_PROBLEM_SEDOV_TAYLOR: ! ------------------------------------- ! ! Subroutine sets the initial conditions for the spherical Sedov-Taylor ! blast problem. ! ! Arguments: ! ! pdata - pointer to the datablock structure of the currently initialized ! block; ! ! References: ! ! [1] Almgren, A. S. et al., ! "CASTRO: A New Compressible Astrophysical Solver. ! I. Hydrodynamics and Self-Gravity", ! The Astrophysical Journal, 2010, vol. 715, pp. 1221-1238, ! http://dx.doi.org/10.1088/0004-637X/715/2/1221 ! !=============================================================================== ! subroutine setup_problem_sedov_taylor(pdata) ! include external procedures and variables ! use blocks , only : block_data use constants , only : d2r #if NDIMS == 3 use constants , only : pi4 #else /* NDIMS == 3 */ use constants , only : pi #endif /* NDIMS == 3 */ use coordinates, only : nn => bcells use coordinates, only : ax, ay, adx, ady, advol #if NDIMS == 3 use coordinates, only : az, adz #endif /* NDIMS == 3 */ use equations , only : prim2cons use equations , only : adiabatic_index use equations , only : nv use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp use parameters , only : get_parameter ! local variables are not implicit by default ! implicit none ! input arguments ! type(block_data), pointer, intent(inout) :: pdata ! default parameter values ! real(kind=8), save :: radius = 1.00d-02 real(kind=8), save :: dens = 1.00d+00 real(kind=8), save :: pres = 1.00d-05 real(kind=8), save :: eexp = 1.00d+00 real(kind=8), save :: buni = 0.00d+00 real(kind=8), save :: angle = 0.00d+00 #if NDIMS == 3 integer , save :: nsubgrid = 10 #endif /* NDIMS == 3 */ ! local saved parameters ! logical , save :: first = .true. real(kind=8), save :: dn_amb, dn_ovr real(kind=8), save :: pr_amb, pr_ovr real(kind=8), save :: bx, by real(kind=8), save :: r2 ! local variables ! integer :: i, j, k #if NDIMS == 3 integer :: ic, jc, kc #endif /* NDIMS == 3 */ real(kind=8) :: xl, yl, xu, yu, rl, ru real(kind=8) :: sn #if NDIMS == 3 real(kind=8) :: zl, zu, dz, dzh real(kind=8) :: xb, yb, zb real(kind=8) :: xt, yt, zt real(kind=8) :: fc_inc #else /* NDIMS == 3 */ real(kind=8) :: rlu, rul real(kind=8) :: xb, yb real(kind=8) :: xt, yt real(kind=8) :: ph #endif /* NDIMS == 3 */ real(kind=8) :: dx, dy, dxh, dyh, dvol real(kind=8) :: fc_amb, fc_ovr ! local arrays ! real(kind=8), dimension(nv,nn) :: q, u real(kind=8), dimension(nn) :: x, y #if NDIMS == 3 real(kind=8), dimension(nn) :: z ! allocatable arrays ! real(kind=8), dimension(:), allocatable :: xm, ym, zm real(kind=8), dimension(:), allocatable :: xp, yp, zp #endif /* NDIMS == 3 */ ! !------------------------------------------------------------------------------- ! ! prepare problem constants during the first subroutine call ! if (first) then ! get problem parameters ! call get_parameter("radius", radius) call get_parameter("dens" , dens ) call get_parameter("pres" , pres ) call get_parameter("eexp" , eexp ) call get_parameter("buni" , buni ) call get_parameter("angle" , angle ) #if NDIMS == 3 ! get the fine grid resolution ! call get_parameter("nsubgrid", nsubgrid) ! correct subgrid resolution if necessary ! nsubgrid = max(1, nsubgrid) #endif /* NDIMS == 3 */ ! calculate the volume of the injection region ! #if NDIMS == 3 dvol = pi4 * radius**3 / 3.0d+00 #else /* NDIMS == 3 */ dvol = pi * radius**2 #endif /* NDIMS == 3 */ ! calculate the overdense and ambient region densities and pressures ! dn_amb = dens dn_ovr = dn_amb pr_amb = pres pr_ovr = (adiabatic_index - 1.0d+00) * eexp / dvol ! calculate initial uniform field components ! if (ibx > 0) then sn = sin(d2r * angle) bx = buni * sqrt(1.0d+00 - sn * sn) by = buni * sn end if ! calculate the square of radius ! r2 = radius * radius ! reset the first execution flag ! first = .false. end if ! first call ! prepare block coordinates ! x(:) = pdata%meta%xmin + ax(pdata%meta%level,:) y(:) = pdata%meta%ymin + ay(pdata%meta%level,:) #if NDIMS == 3 z(:) = pdata%meta%zmin + az(pdata%meta%level,:) #endif /* NDIMS == 3 */ ! calculate mesh intervals and areas ! dx = adx(pdata%meta%level) dy = ady(pdata%meta%level) #if NDIMS == 3 dz = adz(pdata%meta%level) #endif /* NDIMS == 3 */ dxh = 0.5d+00 * dx dyh = 0.5d+00 * dy #if NDIMS == 3 dzh = 0.5d+00 * dz #endif /* NDIMS == 3 */ dvol = advol(pdata%meta%level) #if NDIMS == 3 ! allocate subgrid coordinates ! allocate(xm(nsubgrid), ym(nsubgrid), zm(nsubgrid)) allocate(xp(nsubgrid), yp(nsubgrid), zp(nsubgrid)) ! and generate them ! xm(:) = (1.0d+00 * (/(i, i = 0, nsubgrid - 1)/)) / nsubgrid ym(:) = xm(:) zm(:) = xm(:) xm(:) = xm(:) * dx ym(:) = ym(:) * dy zm(:) = zm(:) * dz xp(:) = (1.0d+00 * (/(i, i = 1, nsubgrid )/)) / nsubgrid yp(:) = xp(:) zp(:) = xp(:) xp(:) = xp(:) * dx yp(:) = yp(:) * dy zp(:) = zp(:) * dz ! calculate the factor increment for the given subgrid ! fc_inc = dvol / nsubgrid**3 #endif /* NDIMS == 3 */ ! set density and pressure of the ambient ! q(idn,:) = dn_amb if (ipr > 0) q(ipr,:) = pr_amb ! reset velocity components ! q(ivx,:) = 0.0d+00 q(ivy,:) = 0.0d+00 q(ivz,:) = 0.0d+00 ! if magnetic field is present, set it to be uniform with the desired strength ! and orientation ! if (ibx > 0) then q(ibx,:) = bx q(iby,:) = by q(ibz,:) = 0.0d+00 q(ibp,:) = 0.0d+00 end if ! iterate over all positions in the YZ plane ! #if NDIMS == 2 k = 1 #endif /* NDIMS == 2 */ #if NDIMS == 3 do k = 1, nn ! calculate the corner Z coordinates ! zl = abs(z(k)) - dzh zu = abs(z(k)) + dzh #endif /* NDIMS == 3 */ do j = 1, nn ! calculate the corner Y coordinates ! yl = abs(y(j)) - dyh yu = abs(y(j)) + dyh ! sweep along the X coordinate ! do i = 1, nn ! calculate the corner X coordinates ! xl = abs(x(i)) - dxh xu = abs(x(i)) + dxh ! calculate the minimum and maximum corner distances from the origin ! #if NDIMS == 3 rl = xl * xl + yl * yl + zl * zl ru = xu * xu + yu * yu + zu * zu #else /* NDIMS == 3 */ rl = xl * xl + yl * yl ru = xu * xu + yu * yu #endif /* NDIMS == 3 */ ! set the initial density and pressure in cells laying completely within ! the blast radius ! if (ru <= r2) then ! set density and pressure for the overpressure region ! q(idn,i) = dn_ovr if (ipr > 0) q(ipr,i) = pr_ovr ! set the initial pressure in the cell completely outside the radius ! else if (rl >= r2) then ! set density and pressure of the ambient ! q(idn,i) = dn_amb if (ipr > 0) q(ipr,i) = pr_amb ! integrate density or pressure in cells which are crossed by the circule with ! the given radius ! else #if NDIMS == 3 ! interpolate the factor using subgrid ! fc_ovr = 0.0d+00 do kc = 1, nsubgrid zb = (zl + zm(kc))**2 zt = (zl + zp(kc))**2 do jc = 1, nsubgrid yb = (yl + ym(jc))**2 yt = (yl + yp(jc))**2 do ic = 1, nsubgrid xb = (xl + xm(ic))**2 xt = (xl + xp(ic))**2 ! update the integration factor depending on the subcell position ! if ((xt + yt + zt) <= r2) then fc_ovr = fc_ovr + fc_inc else if ((xb + yb + zb) < r2) then fc_ovr = fc_ovr + 0.5d+00 * fc_inc end if end do ! ic = 1, nsubgrid end do ! jc = 1, nsubgrid end do ! kc = 1, nsubgrid #else /* NDIMS == 3 */ ! calculate the distance of remaining corners ! rlu = xl * xl + yu * yu rul = xu * xu + yl * yl ! separate in the cases of which corners lay inside, and which outside ! the radius ! if (min(rlu, rul) >= r2) then ! only one cell corner inside the radius ! ! calculate middle coordinates of the radius-edge crossing point ! xb = sqrt(r2 - yl**2) - xl yb = sqrt(r2 - xl**2) - yl ! calculate the sin(½φ), φ, and sin(φ) ! sn = 0.5d+00 * sqrt(xb**2 + yb**2) / radius ph = 2.0d+00 * asin(sn) sn = sin(ph) ! calculate the area of cell intersection with the radius ! fc_ovr = 0.5d+00 * (xb * yb + (ph - sn) * r2) else if (rlu >= r2) then ! two lower corners inside the radius ! ! calculate middle coordinates of the radius-edge crossing point ! yb = sqrt(r2 - xl**2) - yl yt = sqrt(r2 - xu**2) - yl ! calculate the sin(½φ), φ, and sin(φ) ! sn = 0.5d+00 * sqrt(dx**2 + (yt - yb)**2) / radius ph = 2.0d+00 * asin(sn) sn = sin(ph) ! calculate the area of cell intersection with the radius ! fc_ovr = 0.5d+00 * ((yt + yb) * dx + (ph - sn) * r2) else if (rul >= r2) then ! two left corners inside the radius ! ! calculate middle coordinates of the radius-edge crossing point ! xb = sqrt(r2 - yl**2) - xl xt = sqrt(r2 - yu**2) - xl ! calculate the sin(½φ), φ, and sin(φ) ! sn = 0.5d+00 * sqrt((xt - xb)**2 + dy**2) / radius ph = 2.0d+00 * asin(sn) sn = sin(ph) ! calculate the area of cell intersection with the radius ! fc_ovr = 0.5d+00 * ((xt + xb) * dy + (ph - sn) * r2) else ! three corners inside the radius ! ! calculate middle coordinates of the radius-edge crossing point ! xt = xu - sqrt(r2 - yu**2) yt = yu - sqrt(r2 - xu**2) ! calculate the sin(½φ), φ, and sin(φ) ! sn = 0.5d+00 * sqrt(xt**2 + yt**2) / radius ph = 2.0d+00 * asin(sn) sn = sin(ph) ! calculate the area of cell intersection with the radius ! fc_ovr = dvol - 0.5d+00 * (xt * yt - (ph - sn) * r2) end if #endif /* NDIMS == 3 */ ! normalize coefficients ! fc_ovr = fc_ovr / dvol fc_amb = 1.0d+00 - fc_ovr ! integrate density and pressure over the edge cells ! q(idn,i) = fc_ovr * dn_ovr + fc_amb * dn_amb if (ipr > 0) q(ipr,i) = fc_ovr * pr_ovr + fc_amb * pr_amb end if end do ! i = 1, nn ! convert the primitive variables to conservative ones ! call prim2cons(q(:,:), u(:,:)) ! copy the conserved variables to the current block ! pdata%u(:,:,j,k) = u(:,:) ! copy the primitive variables to the current block ! pdata%q(:,:,j,k) = q(:,:) end do ! j = 1, nn #if NDIMS == 3 end do ! k = 1, nn ! deallocate subgrid coordinates ! deallocate(xm, ym, zm) deallocate(xp, yp, zp) #endif /* NDIMS == 3 */ !------------------------------------------------------------------------------- ! end subroutine setup_problem_sedov_taylor ! !=============================================================================== ! ! subroutine SETUP_PROBLEM_IMPLOSION: ! ---------------------------------- ! ! Subroutine sets the initial conditions for the implosion problem. ! ! Arguments: ! ! pdata - pointer to the datablock structure of the currently initialized ! block; ! !=============================================================================== ! subroutine setup_problem_implosion(pdata) ! include external procedures and variables ! use blocks , only : block_data, ndims use constants , only : d2r use coordinates, only : nn => bcells use coordinates, only : ax, ay, adx, ady #if NDIMS == 3 use coordinates, only : az, adz #endif /* NDIMS == 3 */ use equations , only : prim2cons use equations , only : nv use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp use parameters , only : get_parameter ! local variables are not implicit by default ! implicit none ! input arguments ! type(block_data), pointer, intent(inout) :: pdata ! default parameter values ! real(kind=8), save :: sline = 1.50d-01 real(kind=8), save :: adens = 1.00d+00 real(kind=8), save :: apres = 1.00d+00 real(kind=8), save :: drat = 1.25d-01 real(kind=8), save :: prat = 1.40d-01 real(kind=8), save :: buni = 1.00d+00 real(kind=8), save :: bgui = 0.00d+00 real(kind=8), save :: angle = 0.00d+00 ! local saved parameters ! logical , save :: first = .true. real(kind=8), save :: odens = 1.25d-01 real(kind=8), save :: opres = 1.40d-01 ! local variables ! integer :: i, j, k real(kind=8) :: rl, ru, dx, dy, dxh, dyh, ds, dl, dr #if NDIMS == 3 real(kind=8) :: dz, dzh #endif /* NDIMS == 3 */ real(kind=8) :: sn, cs ! local arrays ! real(kind=8), dimension(nv,nn) :: q, u real(kind=8), dimension(nn) :: x, xl, xu real(kind=8), dimension(nn) :: y, yl, yu #if NDIMS == 3 real(kind=8), dimension(nn) :: z, zl, zu #endif /* NDIMS == 3 */ ! !------------------------------------------------------------------------------- ! ! prepare problem constants during the first subroutine call ! if (first) then ! get problem parameters ! call get_parameter("shock_line" , sline ) call get_parameter("ambient_density" , adens ) call get_parameter("ambient_pressure", apres ) call get_parameter("density_ratio" , drat ) call get_parameter("pressure_ratio" , prat ) call get_parameter("buni" , buni ) call get_parameter("bgui" , bgui ) call get_parameter("angle" , angle ) ! calculate parameters ! odens = drat * adens opres = prat * apres ! reset the first execution flag ! first = .false. end if ! first call ! prepare block coordinates ! x(:) = pdata%meta%xmin + ax(pdata%meta%level,:) y(:) = pdata%meta%ymin + ay(pdata%meta%level,:) #if NDIMS == 3 z(:) = pdata%meta%zmin + az(pdata%meta%level,:) #endif /* NDIMS == 3 */ ! calculate mesh intervals and areas ! dx = adx(pdata%meta%level) dy = ady(pdata%meta%level) #if NDIMS == 3 dz = adz(pdata%meta%level) #endif /* NDIMS == 3 */ dxh = 0.5d+00 * dx dyh = 0.5d+00 * dy #if NDIMS == 3 dzh = 0.5d+00 * dz #endif /* NDIMS == 3 */ ! calculate edge coordinates ! xl(:) = abs(x(:)) - dxh xu(:) = abs(x(:)) + dxh yl(:) = abs(y(:)) - dyh yu(:) = abs(y(:)) + dyh #if NDIMS == 3 zl(:) = abs(z(:)) - dzh zu(:) = abs(z(:)) + dzh #endif /* NDIMS == 3 */ ! reset velocity components ! q(ivx,:) = 0.0d+00 q(ivy,:) = 0.0d+00 q(ivz,:) = 0.0d+00 ! if magnetic field is present, set it to be uniform with the desired strength ! and orientation ! if (ibx > 0) then ! calculate the orientation angles ! sn = sin(d2r * angle) cs = sqrt(1.0d+00 - sn * sn) ! set magnetic field components ! q(ibx,:) = buni * cs q(iby,:) = buni * sn q(ibz,:) = bgui q(ibp,:) = 0.0d+00 end if ! iterate over all positions ! #if NDIMS == 2 k = 1 #endif /* NDIMS == 2 */ #if NDIMS == 3 do k = 1, nn #endif /* NDIMS == 3 */ do j = 1, nn do i = 1, nn ! calculate the distance from the origin ! #if NDIMS == 3 rl = xl(i) + yl(j) + zl(k) ru = xu(i) + yu(j) + zu(k) #else /* NDIMS == 3 */ rl = xl(i) + yl(j) ru = xu(i) + yu(j) #endif /* NDIMS == 3 */ ! initialize density and pressure ! if (ru <= sline) then q(idn,i) = odens if (ipr > 0) q(ipr,i) = opres else if (rl >= sline) then q(idn,i) = adens if (ipr > 0) q(ipr,i) = apres else ds = (sline - rl) / dx if (ds <= 1.0d+00) then dl = 5.0d-01 * ds**ndims dr = 1.0d+00 - dl else ds = (ru - sline) / dx dr = 5.0d-01 * ds**ndims dl = 1.0d+00 - dr end if q(idn,i) = adens * dl + odens * dr if (ipr > 0) q(ipr,i) = apres * dl + opres * dr end if end do ! i = 1, im ! convert the primitive variables to conservative ones ! call prim2cons(q(:,:), u(:,:)) ! copy the conserved variables to the current block ! pdata%u(:,:,j,k) = u(:,:) ! copy the primitive variables to the current block ! pdata%q(:,:,j,k) = q(:,:) end do ! j = 1, nn #if NDIMS == 3 end do ! k = 1, nn #endif /* NDIMS == 3 */ !------------------------------------------------------------------------------- ! end subroutine setup_problem_implosion ! !=============================================================================== ! ! subroutine SETUP_PROBLEM_KELVIN_HELMHOLTZ: ! ----------------------------------------- ! ! Subroutine sets the initial conditions for the Kelvin-Helmholtz instability ! problem. ! ! Arguments: ! ! pdata - pointer to the datablock structure of the currently initialized ! block; ! !=============================================================================== ! subroutine setup_problem_kelvin_helmholtz(pdata) ! include external procedures and variables ! use blocks , only : block_data use constants , only : d2r use coordinates, only : nn => bcells use coordinates, only : ay, ady use equations , only : prim2cons use equations , only : nv use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp, isl use parameters , only : get_parameter use random , only : randsym ! local variables are not implicit by default ! implicit none ! input arguments ! type(block_data), pointer, intent(inout) :: pdata ! default parameter values ! real(kind=8), save :: ycut = 2.50d-01 real(kind=8), save :: dens = 1.00d+00 real(kind=8), save :: drat = 2.00d+00 real(kind=8), save :: pres = 2.50d+00 real(kind=8), save :: vamp = 1.00d+00 real(kind=8), save :: vper = 1.00d-02 real(kind=8), save :: buni = 1.00d+00 real(kind=8), save :: bgui = 0.00d+00 real(kind=8), save :: angle = 0.00d+00 ! local saved parameters ! logical , save :: first = .true. ! local variables ! integer :: i, j, k real(kind=8) :: yl, yu, dy, dyh real(kind=8) :: sn, cs ! local arrays ! real(kind=8), dimension(nv,nn) :: q, u real(kind=8), dimension(nn) :: y ! !------------------------------------------------------------------------------- ! ! prepare problem constants during the first subroutine call ! if (first) then ! get problem parameters ! call get_parameter("ycut" , ycut ) call get_parameter("dens" , dens ) call get_parameter("drat" , drat ) call get_parameter("pres" , pres ) call get_parameter("vamp" , vamp ) call get_parameter("vper" , vper ) call get_parameter("buni" , buni ) call get_parameter("bgui" , bgui ) call get_parameter("angle" , angle ) ! reset the first execution flag ! first = .false. end if ! first call ! prepare block coordinates ! y(:) = pdata%meta%ymin + ay(pdata%meta%level,:) ! calculate mesh intervals and areas ! dy = ady(pdata%meta%level) dyh = 0.5d+00 * dy ! set the ambient density and pressure ! q(idn,:) = dens if (ipr > 0) q(ipr,:) = pres ! if magnetic field is present, set it to be uniform with the desired strength ! and orientation ! if (ibx > 0) then ! calculate the orientation angles ! sn = sin(d2r * angle) cs = sqrt(1.0d+00 - sn * sn) ! set magnetic field components ! q(ibx,:) = buni * cs q(iby,:) = buni * sn q(ibz,:) = bgui q(ibp,:) = 0.0d+00 end if ! iterate over all positions in the YZ plane ! #if NDIMS == 2 k = 1 #endif /* NDIMS == 2 */ #if NDIMS == 3 do k = 1, nn #endif /* NDIMS == 3 */ do j = 1, nn ! calculate the corner Y coordinates ! yl = abs(y(j)) - dyh yu = abs(y(j)) + dyh ! set the primitive variables for two regions ! if (yu <= ycut) then q(idn,:) = dens * drat q(ivx,:) = vamp if (isl > 0) q(isl,:) = 1.0d+00 else if (yl >= ycut) then q(idn,:) = dens q(ivx,:) = - vamp if (isl > 0) q(isl,:) = 0.0d+00 else q(idn,:) = dens * ((yu - ycut) * drat + (ycut - yl)) / dy q(ivx,:) = vamp * ((yu - ycut) - (ycut - yl)) / dy if (isl > 0) q(isl,:) = (yu - ycut) / dy end if ! reset remaining velocity components ! q(ivy,:) = 0.0d+00 q(ivz,:) = 0.0d+00 ! set the pressure ! if (ipr > 0) q(ipr,:) = pres ! add a random seed velocity component ! do i = 1, nn q(ivx,i) = q(ivx,i) + vper * randsym() q(ivy,i) = q(ivy,i) + vper * randsym() #if NDIMS == 3 q(ivz,i) = q(ivz,i) + vper * randsym() #endif /* NDIMS == 3 */ end do ! convert the primitive variables to conservative ones ! call prim2cons(q(:,:), u(:,:), .true.) ! copy the conserved variables to the current block ! pdata%u(:,:,j,k) = u(:,:) ! copy the primitive variables to the current block ! pdata%q(:,:,j,k) = q(:,:) end do ! j = 1, nn #if NDIMS == 3 end do ! k = 1, nn #endif /* NDIMS == 3 */ !------------------------------------------------------------------------------- ! end subroutine setup_problem_kelvin_helmholtz ! !=============================================================================== ! ! subroutine SETUP_PROBLEM_RAYLEIGH_TAYLOR: ! ---------------------------------------- ! ! Subroutine sets the initial conditions for the Rayleigh-Taylor instability ! problem. ! ! Arguments: ! ! pdata - pointer to the datablock structure of the currently initialized ! block; ! ! References: ! ! [1] Almgren, A. S. et al., ! "CASTRO: A New Compressible Astrophysical Solver. ! I. Hydrodynamics and Self-Gravity", ! The Astrophysical Journal, 2010, vol. 715, pp. 1221-1238, ! http://dx.doi.org/10.1088/0004-637X/715/2/1221 ! !=============================================================================== ! subroutine setup_problem_rayleigh_taylor(pdata) ! include external procedures and variables ! use blocks , only : block_data use constants , only : pi2, d2r use coordinates, only : xmin, xmax, xlen use coordinates, only : nn => bcells use coordinates, only : ax, ay use equations , only : prim2cons use equations , only : nv use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp, isl use equations , only : csnd2 use parameters , only : get_parameter use random , only : randsym ! local variables are not implicit by default ! implicit none ! input arguments ! type(block_data), pointer, intent(inout) :: pdata ! default parameter values ! real(kind=8), save :: dens = 1.00d+00 real(kind=8), save :: drat = 2.00d+00 real(kind=8), save :: damp = 5.00d-01 real(kind=8), save :: pres = 5.00d+00 real(kind=8), save :: ycut = 0.00d+00 real(kind=8), save :: vper = 0.00d+00 real(kind=8), save :: lper = 1.00d-02 real(kind=8), save :: kper = 1.00d+00 real(kind=8), save :: hdel = 5.00d-03 real(kind=8), save :: gacc = -1.00d+00 real(kind=8), save :: buni = 1.00d+00 real(kind=8), save :: bgui = 0.00d+00 real(kind=8), save :: angle = 0.00d+00 ! local saved parameters ! logical , save :: first = .true. ! local variables ! integer :: i, j, k real(kind=8) :: sn, cs ! local arrays ! real(kind=8), dimension(nv,nn) :: q, u real(kind=8), dimension(nn) :: x, y, yp ! !------------------------------------------------------------------------------- ! ! prepare problem constants during the first subroutine call ! if (first) then ! get problem parameters ! call get_parameter("ycut" , ycut ) call get_parameter("dens" , dens ) call get_parameter("drat" , drat ) call get_parameter("pres" , pres ) call get_parameter("lper" , lper ) call get_parameter("kper" , kper ) call get_parameter("vper" , vper ) call get_parameter("hdel" , hdel ) call get_parameter("gacc" , gacc ) call get_parameter("buni" , buni ) call get_parameter("bgui" , bgui ) call get_parameter("angle" , angle ) ! calculate the density change across the interface ! damp = 5.0d-01 * (drat * dens - dens) ! reset the first execution flag ! first = .false. end if ! first call ! prepare block coordinates ! x(:) = pdata%meta%xmin + ax(pdata%meta%level,:) y(:) = pdata%meta%ymin + ay(pdata%meta%level,:) ! set the ambient density and pressure ! q(idn,:) = dens if (ipr > 0) q(ipr,:) = pres ! if magnetic field is present, set it to be uniform with the desired strength ! and orientation ! if (ibx > 0) then ! calculate the orientation angles ! sn = sin(d2r * angle) cs = sqrt(1.0d+00 - sn * sn) ! set magnetic field components ! q(ibx,:) = buni * cs q(iby,:) = buni * sn q(ibz,:) = bgui q(ibp,:) = 0.0d+00 end if ! prepare density perturbation ! yp(:) = 0.5d+00 * lper * (cos(pi2 * kper * (x(:) - xmin) / xlen) & + cos(pi2 * kper * (xmax - x(:)) / xlen)) + ycut ! iterate over all positions in the YZ plane ! #if NDIMS == 2 k = 1 #endif /* NDIMS == 2 */ #if NDIMS == 3 do k = 1, nn #endif /* NDIMS == 3 */ do j = 1, nn ! set density and pressure ! if (ipr > 0) then if (y(j) <= ycut) then q(idn,:) = dens else q(idn,:) = dens * drat end if q(idn,:) = dens + damp * (1.0d+00 + tanh((y(j) - yp(:)) / hdel)) q(ipr,:) = pres + q(idn,:) * gacc * y(j) else if (y(j) <= ycut) then q(idn,:) = dens * exp(gacc * y(j) / csnd2) else q(idn,:) = dens * drat * exp(gacc * y(j) / csnd2) end if end if if (isl > 0) then if (y(j) <= ycut) then q(isl,:) = 0.0d+00 else q(isl,:) = 1.0d+00 end if end if ! reset the velocity components ! q(ivx,:) = 0.0d+00 q(ivy,:) = 0.0d+00 q(ivz,:) = 0.0d+00 ! add a random seed velocity component ! if (abs(vper) > 0.0d+00) then do i = 1, nn q(ivy,i) = q(ivy,i) + vper * randsym() end do end if ! convert the primitive variables to conservative ones ! call prim2cons(q(:,:), u(:,:), .true.) ! copy the conserved variables to the current block ! pdata%u(:,:,j,k) = u(:,:) ! copy the primitive variables to the current block ! pdata%q(:,:,j,k) = q(:,:) end do ! j = 1, nn #if NDIMS == 3 end do ! k = 1, nn #endif /* NDIMS == 3 */ !------------------------------------------------------------------------------- ! end subroutine setup_problem_rayleigh_taylor ! !=============================================================================== ! ! subroutine SETUP_PROBLEM_ORSZAG_TANG: ! ------------------------------------ ! ! Subroutine sets the initial conditions for the Orszag-Tang test ! problem. ! ! Arguments: ! ! pdata - pointer to the datablock structure of the currently initialized ! block; ! ! References: ! ! [1] Orszag, S. A., Tang, C. -M., ! "Small-scale structure of two-dimensional magnetohydrodynamic ! turbulence", ! Journal of Fluid Mechanics, 1979, 90, pp. 129-143, ! https://doi.org/10.1017/S002211207900210X ! [2] Mignone, A., ! "A simple and accurate Riemann solver for isothermal MHD", ! Journal of Computational Physics, 2007, 225, pp. 1427-1441, ! https://doi.org/10.1016/j.jcp.2007.01.033 ! !=============================================================================== ! subroutine setup_problem_orszag_tang(pdata) use blocks , only : block_data use constants , only : pi2 use coordinates, only : nn => bcells, ax, ay use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp use equations , only : nv, magnetized, csnd, adiabatic_index use equations , only : prim2cons use helpers , only : print_message use parameters , only : get_parameter implicit none type(block_data), pointer, intent(inout) :: pdata real(kind=8), save :: dens = 1.00d+00 real(kind=8), save :: pres = 1.00d+00 real(kind=8), save :: vamp = 1.00d+00 real(kind=8), save :: bamp = 1.00d+00 real(kind=8), save :: bgui = 0.00d+00 logical, save :: first = .true. integer :: j, k real(kind=8), dimension(nv,nn) :: q, u real(kind=8), dimension(nn) :: x, y character(len=*), parameter :: loc = 'PROBLEMS::setup_problem_orszag_tang()' !------------------------------------------------------------------------------- ! if (first) then call get_parameter("dens", dens) call get_parameter("vamp", vamp) call get_parameter("bamp", bamp) call get_parameter("bgui", bgui) if (ipr > 0) then pres = dens / adiabatic_index bamp = bamp / adiabatic_index else vamp = csnd * vamp bamp = csnd * bamp end if first = .false. end if x(:) = pdata%meta%xmin + ax(pdata%meta%level,:) y(:) = pdata%meta%ymin + ay(pdata%meta%level,:) q(idn,:) = dens if (ipr > 0) q(ipr,:) = pres q(ivy,:) = - vamp * sin(pi2 * x(:)) q(ivz,:) = 0.0d+00 if (magnetized) then q(iby,:) = bamp * sin(pi2 * x(:) * 2.0d+00) q(ibz,:) = bgui q(ibp,:) = 0.0d+00 #if NDIMS == 2 k = 1 #endif /* NDIMS == 2 */ #if NDIMS == 3 do k = 1, nn #endif /* NDIMS == 3 */ do j = 1, nn q(ivx,:) = vamp * sin(pi2 * y(j)) q(ibx,:) = bamp * sin(pi2 * y(j)) call prim2cons(q(:,:), u(:,:), .true.) pdata%u(:,:,j,k) = u(:,:) pdata%q(:,:,j,k) = q(:,:) end do #if NDIMS == 3 end do #endif /* NDIMS == 3 */ else call print_message(loc, "Orszag-Tang problem requires MHD!") end if !------------------------------------------------------------------------------- ! end subroutine setup_problem_orszag_tang ! !=============================================================================== ! ! subroutine SETUP_PROBLEM_CURRENT_SHEET: ! -------------------------------------- ! ! Subroutine sets the initial conditions for the current sheet test problem. ! ! Arguments: ! ! pdata - pointer to the datablock structure of the currently initialized ! block; ! !=============================================================================== ! subroutine setup_problem_current_sheet(pdata) ! include external procedures and variables ! use blocks , only : block_data use constants , only : pi2 use coordinates, only : nn => bcells use coordinates, only : ax, ay use equations , only : prim2cons use equations , only : nv, magnetized use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp use parameters , only : get_parameter ! local variables are not implicit by default ! implicit none ! input arguments ! type(block_data), pointer, intent(inout) :: pdata ! default parameter values ! real(kind=8), save :: xpos = 2.50d-01 real(kind=8), save :: xdel = 1.00d-08 real(kind=8), save :: dens = 1.00d+00 real(kind=8), save :: beta = 1.00d-01 real(kind=8), save :: vamp = 1.00d-01 real(kind=8), save :: bamp = 1.00d+00 real(kind=8), save :: bgui = 0.00d+00 ! local saved parameters ! logical, save :: first = .true. ! local variables ! integer :: i, j, k real(kind=8) :: t ! local arrays ! real(kind=8), dimension(nv,nn) :: q, u real(kind=8), dimension(nn) :: x, y, vx ! !------------------------------------------------------------------------------- ! ! retrieve the problem parameters during the first subroutine call ! if (first) then call get_parameter("position" , xpos) call get_parameter("thickness" , xdel) call get_parameter("density" , dens) call get_parameter("beta" , beta) call get_parameter("vel_amplitude", vamp) call get_parameter("mag_amplitude", bamp) call get_parameter("mag_guide" , bgui) first = .false. end if ! prepare the block coordinates ! x(:) = pdata%meta%xmin + ax(pdata%meta%level,:) y(:) = pdata%meta%ymin + ay(pdata%meta%level,:) ! prepare the vector of velocity perturbation (varying along Y) ! vx(:) = vamp * sin(pi2 * y(:)) ! set the fields, along X, which do not vary along the YZ planes ! q(idn,:) = dens q(ivy,:) = 0.0d+00 q(ivz,:) = 0.0d+00 if (ipr > 0) q(ipr,:) = 0.5d+00 * beta ! if magnetic field is present, set its antiparallel configuration ! if (magnetized) then ! set magnetic field components ! q(ibx,:) = 0.0d+00 do i = 1, nn t = (xpos - abs(x(i))) / xdel q(iby,i) = sign(bamp, t) * min(abs(t), 1.0d+00) end do q(ibz,:) = bgui q(ibp,:) = 0.0d+00 ! correct the gas pressure due to the variance of magnetic pressure ! if (ipr > 0) q(ipr,:) = q(ipr,:) + 0.5d+00 * (bamp**2 - q(iby,:)**2) end if ! iterate over all positions in the YZ planes ! #if NDIMS == 2 k = 1 #endif /* NDIMS == 2 */ #if NDIMS == 3 do k = 1, nn #endif /* NDIMS == 3 */ do j = 1, nn ! set the velocity perturbation ! q(ivx,:) = vx(j) ! convert the primitive variables to conservative ones ! call prim2cons(q(:,:), u(:,:)) ! copy the conserved and primitive variables to the current block ! pdata%u(:,:,j,k) = u(:,:) pdata%q(:,:,j,k) = q(:,:) end do ! j = 1, nn #if NDIMS == 3 end do ! k = 1, nn #endif /* NDIMS == 3 */ !------------------------------------------------------------------------------- ! end subroutine setup_problem_current_sheet ! !=============================================================================== ! ! subroutine SETUP_PROBLEM_TEARING: ! -------------------------------- ! ! Subroutine sets the initial conditions for the resistive tearing instability ! test problem. ! ! Arguments: ! ! pdata - pointer to the datablock structure of the currently initialized ! block; ! ! References: ! ! [1] Lapenta, G., ! "Self-Feeding Turbulent Magnetic Reconnection on Macroscopic Scales", ! Physical Review Letters, 2008, vol. 100, id. 235001, ! http://dx.doi.org/10.1103/PhysRevLett.100.235001 ! ! [2] Landi, S., Del Zanna, L., Papini, E., Pucci, F., and Velli, M. ! "Resistive Magnetohydrodynamics Simulations of the Ideal Tearing Mode", ! The Astrophysical Journal, 2015, vol. 806, id. 131, ! http://dx.doi.org/10.1088/0004-637X/806/1/131 ! !=============================================================================== ! subroutine setup_problem_tearing(pdata) ! include external procedures and variables ! use blocks , only : block_data use constants , only : pi2 use coordinates, only : nn => bcells use coordinates, only : ax, ay, adx, ady use equations , only : prim2cons use equations , only : nv, magnetized, csnd2 use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp use parameters , only : get_parameter ! local variables are not implicit by default ! implicit none ! input arguments ! type(block_data), pointer, intent(inout) :: pdata ! default parameter values ! real(kind=8), save :: beta = 1.00d-01 real(kind=8), save :: delta = 1.00d-08 real(kind=8), save :: eta = 1.00d-08 real(kind=8), save :: zeta = 0.00d+00 real(kind=8), save :: dens = 1.00d+00 real(kind=8), save :: pres = 5.00d-02 real(kind=8), save :: pmag = 1.00d+00 real(kind=8), save :: bamp = 1.00d+00 real(kind=8), save :: bnor = 0.00d+00 real(kind=8), save :: bgui = 0.00d+00 real(kind=8), save :: bper = 1.00d-02 real(kind=8), save :: vper = 1.00d-03 real(kind=8), save :: kper = 1.00d+00 real(kind=8), save :: ss = 1.00d+00 ! local saved parameters ! logical, save :: first = .true. ! local variables ! integer :: i, j, k real(kind=8) :: dx, dy, t ! local arrays ! real(kind=8), dimension(nv,nn) :: q, u real(kind=8), dimension(nn) :: xc, xl, xu, yl, yu, xi real(kind=8), dimension(nn) :: vx, vy, bx, by, b0 ! !------------------------------------------------------------------------------- ! ! retrieve the problem parameters during the first subroutine call ! if (first) then call get_parameter("resistivity", eta) call get_parameter("beta" , beta) call get_parameter("zeta" , zeta) call get_parameter("density" , dens) call get_parameter("bamp" , bamp) call get_parameter("bnor" , bnor) call get_parameter("bgui" , bgui) call get_parameter("bper" , bper) call get_parameter("vper" , vper) call get_parameter("kper" , kper) ss = abs(bamp) / (sqrt(dens) * max(1.0d-16, eta)) delta = ss**(-1.0d+00/3.0d+00) if (abs(bper) > 0.0d+00) vper = 0.0d+00 pres = 0.5d+00 * beta pmag = 0.5d+00 * (bamp**2 + bnor**2 + bgui**2) first = .false. end if ! prepare the block coordinates ! dx = adx(pdata%meta%level) dy = ady(pdata%meta%level) xc(:) = pdata%meta%xmin + ax(pdata%meta%level,:) xl(:) = pdata%meta%xmin + ax(pdata%meta%level,:) - 0.5d+00 * dx xu(:) = xl(:) + dx yl(:) = pdata%meta%ymin + ay(pdata%meta%level,:) - 0.5d+00 * dy yu(:) = yl(:) + dy ! prepare the vector of velocity perturbation (varying along Y) ! if (abs(vper) > 0.0d+00) then xi(:) = xc(:) * sqrt(ss) vx(:) = vper * tanh(xi(:)) * exp(- xi(:)**2) vy(:) = vper * (2.0d+00 * xi(:) * tanh(xi(:)) - 1.0d+00 / cosh(xi(:))**2)& * exp(- xi(:)**2) * sqrt(ss) / (pi2 * kper) end if if (abs(bper) > 0.0d+00) then bx(:) = bper * (sin(pi2 * xl(:)) - sin(pi2 * xu(:))) / (pi2 * dx) by(:) = bper * (cos(pi2 * xl(:)) - cos(pi2 * xu(:))) / (pi2 * dx * kper) end if ! set the fields, along X, which do not vary along the YZ planes ! q(idn,:) = dens q(ivx,:) = 0.0d+00 q(ivy,:) = 0.0d+00 q(ivz,:) = 0.0d+00 if (ipr > 0) q(ipr,:) = pres ! if magnetic field is present, set its antiparallel configuration ! if (magnetized) then ! set magnetic field components ! q(ibx,:) = bnor do i = 1, nn t = delta * (log(cosh(xu(i) / delta)) & - log(cosh(xl(i) / delta))) / dx q(iby,i) = bamp * max(-1.0d+00, min(1.0d+00, t)) q(ibz,i) = zeta * sqrt(bamp**2 - q(iby,i)**2) + bgui end do q(ibp,:) = 0.0d+00 b0(:) = q(iby,:) ! correct the gas pressure due to the variance of magnetic pressure ! if (ipr > 0) then q(ipr,:) = pres + pmag - 0.5d+00 * sum(q(ibx:ibz,:)**2,1) else q(idn,:) = (pres + pmag - 0.5d+00 * sum(q(ibx:ibz,:)**2,1)) / csnd2 end if end if ! iterate over all positions in the YZ planes ! #if NDIMS == 2 k = 1 #endif /* NDIMS == 2 */ #if NDIMS == 3 do k = 1, nn #endif /* NDIMS == 3 */ do j = 1, nn ! set the perturbation ! if (abs(vper) > 0.0d+00) then q(ivx,:) = vx(:) * (sin(pi2 * kper * yl(j)) & - sin(pi2 * kper * yu(j))) / (pi2 * kper * dy) q(ivy,:) = vy(:) * (cos(pi2 * kper * yu(j)) & - cos(pi2 * kper * yl(j))) / (pi2 * kper * dy) end if if (abs(bper) > 0.0d+00) then q(ibx,:) = bx(:) * (cos(pi2 * kper * yu(j)) & - cos(pi2 * kper * yl(j))) / (pi2 * kper * dy) & + bnor q(iby,:) = by(:) * (sin(pi2 * kper * yl(j)) & - sin(pi2 * kper * yu(j))) / (pi2 * kper * dy) & + b0(:) end if ! convert the primitive variables to conservative ones ! call prim2cons(q(:,:), u(:,:)) ! copy the conserved and primitive variables to the current block ! pdata%u(:,:,j,k) = u(:,:) pdata%q(:,:,j,k) = q(:,:) end do ! j = 1, nn #if NDIMS == 3 end do ! k = 1, nn #endif /* NDIMS == 3 */ !------------------------------------------------------------------------------- ! end subroutine setup_problem_tearing ! !=============================================================================== ! ! subroutine SETUP_PROBLEM_TURBULENCE: ! ----------------------------------- ! ! Subroutine sets the initial conditions for the turbulence problem. ! ! Arguments: ! ! pdata - pointer to the datablock structure of the currently initialized ! block; ! !=============================================================================== ! subroutine setup_problem_turbulence(pdata) use blocks , only : block_data use coordinates, only : nn => bcells use equations , only : nv, magnetized, adiabatic_index, csnd2 use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp use equations , only : prim2cons use parameters , only : get_parameter implicit none type(block_data), pointer, intent(inout) :: pdata real(kind=8), save :: beta = 1.00d+00 real(kind=8), save :: dens = 1.00d+00 real(kind=8), save :: bamp = 1.00d+00 real(kind=8), save :: bgui = 0.00d+00 real(kind=8), save :: pres = 5.00d-02 real(kind=8), save :: pmag = 1.00d+00 logical, save :: first = .true. integer :: j, k real(kind=8), dimension(nv,nn) :: q, u !------------------------------------------------------------------------------- ! if (first) then call get_parameter("density", dens) pres = csnd2 * dens if (ipr > 0) pres = pres / adiabatic_index if (magnetized) then call get_parameter("bamp", bamp) call get_parameter("bgui", bgui) pmag = 0.5d+00 * (bamp**2 + bgui**2) beta = pres / pmag end if first = .false. end if q(idn,:) = dens q(ivx,:) = 0.0d+00 q(ivy,:) = 0.0d+00 q(ivz,:) = 0.0d+00 if (ipr > 0) q(ipr,:) = pres if (magnetized) then q(ibx,:) = bamp q(iby,:) = 0.0d+00 q(ibz,:) = bgui q(ibp,:) = 0.0d+00 end if call prim2cons(q(:,:), u(:,:)) #if NDIMS == 2 k = 1 #endif /* NDIMS == 2 */ #if NDIMS == 3 do k = 1, nn #endif /* NDIMS == 3 */ do j = 1, nn pdata%u(:,:,j,k) = u(:,:) pdata%q(:,:,j,k) = q(:,:) end do #if NDIMS == 3 end do #endif /* NDIMS == 3 */ !------------------------------------------------------------------------------- ! end subroutine setup_problem_turbulence ! !=============================================================================== ! ! subroutine SETUP_PROBLEM_WIND: ! ----------------------------- ! ! Subroutine sets the initial conditions for the stellar wind problem. ! ! Arguments: ! ! pdata - pointer to the datablock structure of the currently initialized ! block; ! !=============================================================================== ! subroutine setup_problem_wind(pdata) use blocks , only : block_data use coordinates, only : nn => bcells use coordinates, only : ax, ay #if NDIMS == 3 use coordinates, only : az #endif /* NDIMS == 3 */ use equations , only : nv, magnetized, adiabatic_index, csnd, csnd2 use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp use equations , only : prim2cons use parameters , only : get_parameter implicit none type(block_data), pointer, intent(inout) :: pdata real(kind=8), save :: radius = 1.00d-01 real(kind=8), save :: dens = 1.00d+06 real(kind=8), save :: sonic = 3.00d+00 real(kind=8), save :: wind = 5.00d-01 real(kind=8), save :: pres = 5.00d-02 real(kind=8), save :: bamp = 1.00d-03 logical, save :: first = .true. integer :: i, j, k real(kind=8) :: r, r2, s, w real(kind=8), dimension(nv,nn) :: q, u real(kind=8), dimension(nn) :: x, y #if NDIMS == 3 real(kind=8), dimension(nn) :: z #endif /* NDIMS == 3 */ !------------------------------------------------------------------------------- ! if (first) then call get_parameter("radius" , radius) call get_parameter("density" , dens ) call get_parameter("sonic" , sonic ) call get_parameter("wind_speed", wind ) if (ipr > 0) then pres = (wind / sonic)**2 * dens / adiabatic_index else csnd = wind / sonic csnd2 = csnd * csnd end if if (magnetized) & call get_parameter("bamp", bamp) first = .false. end if x(:) = pdata%meta%xmin + ax(pdata%meta%level,:) y(:) = pdata%meta%ymin + ay(pdata%meta%level,:) #if NDIMS == 3 z(:) = pdata%meta%zmin + az(pdata%meta%level,:) #endif /* NDIMS == 3 */ #if NDIMS == 3 do k = 1, nn #else /* NDIMS == 3 */ k = 1 #endif /* NDIMS == 3 */ do j = 1, nn do i = 1, nn #if NDIMS == 3 r2 = x(i)**2 + y(j)**2 + z(k)**2 #else /* NDIMS == 3 */ r2 = x(i)**2 + y(j)**2 #endif /* NDIMS == 3 */ r = max(sqrt(r2), 1.0d-16) s = max(1.0d+00, r / radius)**(NDIMS - 1) q(idn,i) = dens / s if (ipr > 0) q(ipr,i) = pres / s**adiabatic_index q(ivx,i) = wind * x(i) / r q(ivy,i) = wind * y(j) / r #if NDIMS == 3 q(ivz,i) = wind * z(k) / r #else /* NDIMS == 3 */ q(ivz,i) = 0.0d+00 #endif /* NDIMS == 3 */ if (magnetized) then w = max(1.0d+00, r / radius)**(5.0d-1 * NDIMS) q(ibx,i) = 0.0d+00 q(iby,i) = 0.0d+00 q(ibz,i) = bamp / w q(ibp,i) = 0.0d+00 end if end do call prim2cons(q(:,:), u(:,:)) pdata%q(:,:,j,k) = q(:,:) pdata%u(:,:,j,k) = u(:,:) end do #if NDIMS == 3 end do #endif /* NDIMS == 3 */ !------------------------------------------------------------------------------- ! end subroutine setup_problem_wind !=============================================================================== ! end module problems