!!****************************************************************************** !! !! 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: SHAPES !! !! This modules handles the update of the shapes embedded in the domain. In !! such a region, the primitive and conservative variables have to be updated !! after each temporal integration substep. !! !!****************************************************************************** ! module shapes implicit none ! interfaces for procedure to update custom shapes ! abstract interface subroutine update_shapes_iface(pdata, time, dt) use blocks, only : block_data implicit none type(block_data), pointer, intent(inout) :: pdata real(kind=8) , intent(in) :: time, dt end subroutine end interface ! pointer to the shape update subroutine ! procedure(update_shapes_iface), pointer, save :: update_shapes => null() ! a flag to indicate if shapes are on ! logical, save :: enabled = .false. private public :: initialize_shapes, finalize_shapes, print_shapes public :: update_shapes !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! contains ! !=============================================================================== !! !!*** PUBLIC SUBROUTINES ***************************************************** !! !=============================================================================== ! !=============================================================================== ! ! subroutine INITIALIZE_SHAPES: ! ---------------------------- ! ! Subroutine initializes module SHAPES. ! ! Arguments: ! ! status - an integer flag for error return value; ! !=============================================================================== ! subroutine initialize_shapes(status) use parameters, only : get_parameter implicit none integer, intent(out) :: status character(len=64) :: problem_name = "blast" character(len=64) :: enable_shapes = "off" !------------------------------------------------------------------------------- ! status = 0 ! get the problem name ! call get_parameter("problem" , problem_name ) call get_parameter("enable_shapes", enable_shapes) ! set the correct procedure pointer if shapes are enabled ! select case(trim(enable_shapes)) case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES") ! turn on the enable shape flag ! enabled = .true. ! select the shape update subroutine depending on the problem ! select case(trim(problem_name)) ! general test problems ! case("blast") update_shapes => update_shapes_blast case("stellar_wind", "stellar-wind", "wind") update_shapes => update_shapes_wind end select case default ! by default the shape update is turned off, so reset the procedure pointer ! update_shapes => update_shapes_none end select !------------------------------------------------------------------------------- ! end subroutine initialize_shapes ! !=============================================================================== ! ! subroutine FINALIZE_SHAPES: ! -------------------------- ! ! Subroutine releases memory used by the module. ! ! Arguments: ! ! status - an integer flag for error return value; ! !=============================================================================== ! subroutine finalize_shapes(status) implicit none integer, intent(out) :: status !------------------------------------------------------------------------------- ! status = 0 nullify(update_shapes) !------------------------------------------------------------------------------- ! end subroutine finalize_shapes ! !=============================================================================== ! ! subroutine PRINT_SHAPES: ! ----------------------- ! ! Subroutine prints module parameters and settings. ! ! Arguments: ! ! verbose - a logical flag turning the information printing; ! !=============================================================================== ! subroutine print_shapes(verbose) use helpers, only : print_section, print_parameter implicit none logical, intent(in) :: verbose !------------------------------------------------------------------------------- ! if (verbose) call print_parameter(verbose, "embedded shapes", enabled, "on") !------------------------------------------------------------------------------- ! end subroutine print_shapes ! !=============================================================================== !! !!*** PRIVATE SUBROUTINES **************************************************** !! !=============================================================================== ! !=============================================================================== ! ! subroutine UPDATE_SHAPES_NONE: ! ----------------------------- ! ! Subroutine does not do anything, but it is required to define the interface ! for other specific shape update subroutines. ! ! Arguments: ! ! pdata - pointer to the data block structure of the currently initialized ! block; ! time - time at the moment of update; ! dt - time step since the last update; ! !=============================================================================== ! subroutine update_shapes_none(pdata, time, dt) use blocks, only : block_data implicit none type(block_data), pointer, intent(inout) :: pdata real(kind=8) , intent(in) :: time, dt !------------------------------------------------------------------------------- ! !------------------------------------------------------------------------------- ! end subroutine update_shapes_none ! !=============================================================================== ! ! subroutine UPDATE_SHAPES_BLAST: ! ------------------------------ ! ! Subroutine resets the primitive and conserved variables within a defined ! shape for the blast problem. ! ! Arguments: ! ! pdata - pointer to the data block structure of the currently initialized ! block; ! time - time at the moment of update; ! dt - time step since the last update; ! !=============================================================================== ! subroutine update_shapes_blast(pdata, time, dt) use blocks , only : block_data 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 : adiabatic_index use equations , only : nv use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp use parameters , only : get_parameter implicit none type(block_data), pointer, intent(inout) :: pdata real(kind=8) , intent(in) :: time, dt 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 logical , save :: first = .true. !$omp threadprivate(first) real(kind=8), save :: r2 real(kind=8), save :: dn_ovr, pr_ovr, bx_ovr, by_ovr integer :: i, j, k real(kind=8) :: xl, yl, xu, yu, rl, ru real(kind=8) :: dx, dy, dxh, dyh, daxy #if NDIMS == 3 real(kind=8) :: zl, zu, dz, dzh #else /* NDIMS == 3 */ real(kind=8) :: xb, yb, xt, yt real(kind=8) :: fc_amb, fc_ovr #endif /* NDIMS == 3 */ real(kind=8), dimension(nv) :: qi 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 NDIMS == 2 k = 1 #endif /* NDIMS == 2 */ 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 ) ! set the conditions inside the radius ! if (ipr > 0) then dn_ovr = dens pr_ovr = dens * ratio * csnd * csnd / adiabatic_index else dn_ovr = dens * ratio end if bx_ovr = buni * cos(d2r * angle) by_ovr = buni * sin(d2r * angle) ! 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 */ daxy = dx * dy ! set the conditions inside the radius ! if (ipr > 0) then qi(idn) = dn_ovr qi(ipr) = pr_ovr else qi(idn) = dn_ovr end if qi(ivx) = 0.0d+00 qi(ivy) = 0.0d+00 qi(ivz) = 0.0d+00 if (ibx > 0) then qi(ibx) = bx_ovr qi(iby) = by_ovr qi(ibz) = 0.0d+00 qi(ibp) = 0.0d+00 end if ! iterate over all positions in the YZ plane ! #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 primitive variables ! q(1:nv,i) = qi(1:nv) ! variables in cells completely outside the given radius are not changed ! else if (rl >= r2) then q(1:nv,i) = pdata%q(1:nv,i,j,k) ! integrate density or pressure in cells which are crossed by the circle with ! the given radius ! else #if NDIMS == 3 ! in 3D simply set the ambient values since the integration is more complex ! q(1:nv,i) = pdata%q(1:nv,i,j,k) #else /* NDIMS == 3 */ ! calculate the bounds of area integration ! xb = max(xl, sqrt(max(0.0d+00, r2 - yu * yu))) xt = min(xu, sqrt(max(0.0d+00, r2 - yl * yl))) yb = max(yl, sqrt(max(0.0d+00, r2 - xu * xu))) yt = min(yu, sqrt(max(0.0d+00, r2 - xl * xl))) ! integrate the area below the circle within the current cell for both ! functions, y = f(x) and x = g(y), and then average them to be sure that we ! are getting the ideal symmetry ! fc_ovr = 0.5d+00 * (r2 * (asin(xt / radius) - asin(xb / radius)) & + (xt * yb - xb * yt)) - yl * (xt - xb) fc_ovr = fc_ovr + (xb - xl) * dy fc_amb = 0.5d+00 * (r2 * (asin(yt / radius) - asin(yb / radius)) & + (yt * xb - yb * xt)) - xl * (yt - yb) fc_amb = fc_amb + (yb - yl) * dx fc_ovr = 0.5d+00 * (fc_ovr + fc_amb) ! normalize coefficients ! fc_ovr = fc_ovr / daxy fc_amb = 1.0d+00 - fc_ovr ! integrate the primitive variables over the edge cells ! q(1:nv,i) = fc_ovr * qi(1:nv) + fc_amb * pdata%q(1:nv,i,j,k) #endif /* NDIMS == 3 */ 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 #endif /* NDIMS == 3 */ !------------------------------------------------------------------------------- ! end subroutine update_shapes_blast ! !=============================================================================== ! ! subroutine UPDATE_SHAPES_WIND: ! ----------------------------- ! ! Subroutine resets the primitive and conserved variables within a defined ! shape for the wind problem. ! ! Arguments: ! ! pdata - pointer to the data block structure of the currently initialized ! block; ! time - time at the moment of update; ! dt - time step since the last update; ! !=============================================================================== ! subroutine update_shapes_wind(pdata, time, dt) 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 : prim2cons use equations , only : adiabatic_index, csnd, csnd2 use equations , only : nv, magnetized use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp use parameters , only : get_parameter implicit none type(block_data), pointer, intent(inout) :: pdata real(kind=8) , intent(in) :: time, dt 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. !$omp threadprivate(first) real(kind=8), save :: r2, r integer :: i, j, k 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("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) if (r <= radius) then pdata%q(idn,i,j,k) = dens if (ipr > 0) pdata%q(ipr,i,j,k) = pres pdata%q(ivx,i,j,k) = wind * x(i) / r pdata%q(ivy,i,j,k) = wind * y(j) / r #if NDIMS == 3 pdata%q(ivz,i,j,k) = wind * z(k) / r #else /* NDIMS == 3 */ pdata%q(ivz,i,j,k) = 0.0d+00 #endif /* NDIMS == 3 */ if (magnetized) then pdata%q(ibx,i,j,k) = 0.0d+00 pdata%q(iby,i,j,k) = 0.0d+00 pdata%q(ibz,i,j,k) = bamp pdata%q(ibp,i,j,k) = 0.0d+00 end if call prim2cons(pdata%q(:,i,j,k), pdata%u(:,i,j,k)) end if end do end do #if NDIMS == 3 end do #endif /* NDIMS == 3 */ !------------------------------------------------------------------------------- ! end subroutine update_shapes_wind !=============================================================================== ! end module shapes