gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'master' into reconnection

Browse Source 
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2021-11-26 17:01:21 -03:00
commit 93e5ddbea4
7 changed files with 163 additions and 138 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

155
sources/boundaries.F90
View File

@ -35,6 +35,18 @@ module boundaries
implicit none
! interfaces for custom boundary conditions
!
abstract interface
subroutine custom_boundary_iface(idir, il, iu, jl, ju, t, dt, x, y, z, qn)
implicit none
integer , intent(in) :: idir, il, iu, jl, ju
real(kind=8) , intent(in) :: t, dt
real(kind=8), dimension(:) , intent(in) :: x, y, z
real(kind=8), dimension(:,:,:,:), intent(inout) :: qn
end subroutine
end interface
! parameters corresponding to the boundary type
!
integer, parameter :: bnd_periodic = 0
@ -44,6 +56,13 @@ module boundaries
integer, parameter :: bnd_gravity = 4
integer, parameter :: bnd_user = 5
! pointers to the custom boundary conditions, one per direction
!
procedure(custom_boundary_iface), pointer, save :: &
custom_boundary_x => null(), &
custom_boundary_y => null(), &
custom_boundary_z => null()
! variable to store boundary type flags
!
integer, dimension(3,2), save :: bnd_type = bnd_periodic
@ -56,14 +75,11 @@ module boundaries
integer , dimension(:,:), allocatable, save :: bcount
#endif /* MPI */
! by default everything is private
!
private
! declare public subroutines
!
public :: initialize_boundaries, finalize_boundaries, print_boundaries
public :: boundary_variables, boundary_fluxes
public :: custom_boundary_x,custom_boundary_y, custom_boundary_z
public :: bnd_type, bnd_periodic
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
@ -357,28 +373,22 @@ module boundaries
! Arguments:
!
! t, dt - time and time increment;
! status - the call status;
!
!===============================================================================
!
subroutine boundary_variables(t, dt)
subroutine boundary_variables(t, dt, status)
! import external procedures and variables
!
use blocks , only : ndims
use coordinates, only : minlev, maxlev
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), intent(in) :: t, dt
integer , intent(out) :: status
! local variables
!
integer :: idir
!
!-------------------------------------------------------------------------------
!
#if NDIMS == 3
@ -423,7 +433,7 @@ module boundaries
! update specific boundaries
!
call boundaries_specific(t, dt)
call boundaries_specific(t, dt, status)
#if NDIMS == 3
! prolong face boundaries from lower level blocks
@ -447,7 +457,7 @@ module boundaries
! update specific boundaries
!
call boundaries_specific(t, dt)
call boundaries_specific(t, dt, status)
! convert updated primitive variables to conservative ones in all ghost cells
!
@ -1370,13 +1380,12 @@ module boundaries
! Arguments:
!
! t, dt - time and time increment;
! status - the call status;
!
!===============================================================================
!
subroutine boundaries_specific(t, dt)
subroutine boundaries_specific(t, dt, status)
! import external procedures and variables
!
use blocks , only : block_meta, block_leaf
use blocks , only : list_leaf
use blocks , only : ndims, nsides
@ -1390,52 +1399,32 @@ module boundaries
use mpitools , only : nproc
#endif /* MPI */
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), intent(in) :: t, dt
integer , intent(out) :: status
! local pointers
!
type(block_meta), pointer :: pmeta
type(block_leaf), pointer :: pleaf
! local variables
!
integer :: i, j, k = 1, n
#if NDIMS == 2
integer :: m
#endif /* NDIMS == 2 */
! local arrays
!
real(kind=8), dimension(nn) :: x
real(kind=8), dimension(nn) :: y
real(kind=8), dimension(nn) :: x, y
#if NDIMS == 3
real(kind=8), dimension(nn) :: z
#else /* NDIMS == 3 */
real(kind=8), dimension( 1) :: z
real(kind=8), dimension( 1) :: z = 0.0d+00
#endif /* NDIMS == 3 */
!
!-------------------------------------------------------------------------------
!
! associate pleaf with the first block on the leaf list
!
pleaf => list_leaf
! scan all leaf meta blocks in the list
!
do while(associated(pleaf))
! get the associated meta block
!
pmeta => pleaf%meta
! process only if this block is marked for update
!
if (pmeta%update) then
#ifdef MPI
@ -1450,8 +1439,6 @@ module boundaries
y(:) = pmeta%ymin + ay(pmeta%level,:)
#if NDIMS == 3
z(:) = pmeta%zmin + az(pmeta%level,:)
#else /* NDIMS == 3 */
z( : ) = 0.0d+00
#endif /* NDIMS == 3 */
#if NDIMS == 2
@ -1478,7 +1465,7 @@ module boundaries
if (.not. associated(pmeta%edges(i,j,m)%ptr)) &
call block_boundary_specific(n, (/ i, j, k /) &
, t, dt, x(:), y(:), z(:) &
, pmeta%data%q(:,:,:,:))
, pmeta%data%q(:,:,:,:), status)
end do ! i = 1, sides
end do ! j = 1, sides
@ -1507,7 +1494,7 @@ module boundaries
if (.not. associated(pmeta%faces(i,j,k,n)%ptr)) &
call block_boundary_specific(n, (/ i, j, k /) &
, t, dt, x(:), y(:), z(:) &
, pmeta%data%q(:,:,:,:))
, pmeta%data%q(:,:,:,:), status)
end do ! i = 1, sides
end do ! j = 1, sides
@ -1524,11 +1511,8 @@ module boundaries
end if ! if pmeta marked for update
! associate pleaf with the next leaf on the list
!
pleaf => pleaf%next
end do ! over leaf blocks
end do
!-------------------------------------------------------------------------------
!
@ -5081,10 +5065,11 @@ module boundaries
! t, dt - time and time increment;
! x, y, z - the block coordinates;
! qn - the variable array;
! status - the call status;
!
!===============================================================================
!
subroutine block_boundary_specific(nc, side, t, dt, x, y, z, qn)
subroutine block_boundary_specific(nc, side, t, dt, x, y, z, qn, status)
use coordinates , only : nn => bcells, nh => bcells_half, ng => nghosts
use coordinates , only : nb, ne, nbl, neu
@ -5096,18 +5081,15 @@ module boundaries
use equations , only : csnd2
use gravity , only : gravitational_acceleration
use helpers , only : print_message
use user_problem, only : user_boundary_x, user_boundary_y, &
user_boundary_z
implicit none
integer , intent(in) :: nc
integer , dimension(3) , intent(in) :: side
real(kind=8) , intent(in) :: t, dt
real(kind=8), dimension(:) , intent(inout) :: x
real(kind=8), dimension(:) , intent(inout) :: y
real(kind=8), dimension(:) , intent(inout) :: z
real(kind=8), dimension(:) , intent(inout) :: x, y, z
real(kind=8), dimension(:,:,:,:), intent(inout) :: qn
integer , intent(out) :: status
integer :: i, il, iu, is, it, im1, ip1
integer :: j, jl, ju, js, jt, jm1, jp1
@ -5124,6 +5106,8 @@ module boundaries
!-------------------------------------------------------------------------------
!
status = 0
! apply specific boundaries depending on the direction
!
select case(nc)
@ -5288,18 +5272,25 @@ module boundaries
!
case(bnd_user)
call user_boundary_x(side(1), jl, ju, kl, ku &
, t, dt, x(:), y(:), z(:), qn(:,:,:,:))
if (associated(custom_boundary_z)) then
call custom_boundary_x(side(1), jl, ju, kl, ku, &
t, dt, x(:), y(:), z(:), qn(:,:,:,:))
else
call print_message(loc, &
"No subroutine associated with custom_boundary_x()!")
status = 1
return
end if
! wrong boundary conditions
!
case default
if (side(1) == 1) then
call print_message(loc, "Wrong left X boundary type!")
call print_message(loc, "Wrong the left X-boundary type!")
else
call print_message(loc, "Wrong right X boundary type!")
call print_message(loc, "Wrong the right X-boundary type!")
end if
status = 1
return
end select
@ -5464,18 +5455,25 @@ module boundaries
!
case(bnd_user)
call user_boundary_y(side(2), il, iu, kl, ku &
, t, dt, x(:), y(:), z(:), qn(:,:,:,:))
if (associated(custom_boundary_z)) then
call custom_boundary_y(side(2), il, iu, kl, ku, &
t, dt, x(:), y(:), z(:), qn(:,:,:,:))
else
call print_message(loc, &
"No subroutine associated with custom_boundary_y()!")
status = 1
return
end if
! wrong boundary conditions
!
case default
if (side(2) == 1) then
call print_message(loc, "Wrong left Y boundary type!")
call print_message(loc, "Wrong the left Y-boundary type!")
else
call print_message(loc, "Wrong right Y boundary type!")
call print_message(loc, "Wrong the right Y-boundary type!")
end if
status = 1
return
end select
@ -5636,18 +5634,25 @@ module boundaries
!
case(bnd_user)
call user_boundary_z(side(3), il, iu, jl, ju &
, t, dt, x(:), y(:), z(:), qn(:,:,:,:))
if (associated(custom_boundary_z)) then
call custom_boundary_z(side(3), il, iu, jl, ju, &
t, dt, x(:), y(:), z(:), qn(:,:,:,:))
else
call print_message(loc, &
"No subroutine associated with custom_boundary_z()!")
status = 1
return
end if
! wrong boundary conditions
!
case default
if (side(3) == 1) then
call print_message(loc, "Wrong left Z boundary type!")
call print_message(loc, "Wrong the left Z-boundary type!")
else
call print_message(loc, "Wrong right Z boundary type!")
call print_message(loc, "Wrong the right Z-boundary type!")
end if
status = 1
return
end select

5
sources/evolution.F90
View File

@ -1026,7 +1026,6 @@ module evolution
!
call print_sources(verbose)
call print_forcing(verbose)
call print_shapes(verbose)
if (verbose) then
@ -1050,6 +1049,7 @@ module evolution
end if
call print_shapes(verbose)
call print_schemes(verbose)
!-------------------------------------------------------------------------------
@ -3751,7 +3751,8 @@ module evolution
#endif /* MPI */
if (status /= 0) go to 200
call boundary_variables(tm, dtm)
call boundary_variables(tm, dtm, status)
if (status /= 0) go to 200
pdata => list_data
do while (associated(pdata))

8
sources/problems.F90
View File

@ -70,7 +70,6 @@ module problems
use mesh , only : setup_problem
use parameters , only : get_parameter
use user_problem, only : initialize_user_problem, setup_user_problem
implicit none
@ -112,9 +111,6 @@ module problems
setup_problem => setup_problem_tearing
case default
setup_problem => setup_user_problem
call initialize_user_problem(problem, rcount, verbose, status)
end select
@ -137,8 +133,6 @@ module problems
!
subroutine finalize_problems(status)
use user_problem, only : finalize_user_problem
implicit none
integer, intent(out) :: status
@ -147,8 +141,6 @@ module problems
!
status = 0
call finalize_user_problem(status)
!-------------------------------------------------------------------------------
!
end subroutine finalize_problems

25
sources/shapes.F90
View File

@ -33,11 +33,12 @@ module shapes
implicit none
! interfaces for procedure pointers
! 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
@ -51,12 +52,8 @@ module shapes
!
logical, save :: enabled = .false.
! by default everything is private
!
private
! declare public subroutines
!
public :: initialize_shapes, finalize_shapes, print_shapes
public :: update_shapes
@ -85,8 +82,7 @@ module shapes
!
subroutine initialize_shapes(status)
use parameters , only : get_parameter
use user_problem, only : update_user_shapes
use parameters, only : get_parameter
implicit none
@ -122,11 +118,6 @@ module shapes
case("blast")
update_shapes => update_shapes_blast
! no shape update
!
case default
update_shapes => update_user_shapes
end select
case default
@ -164,8 +155,6 @@ module shapes
!
status = 0
! nullify procedure pointers
!
nullify(update_shapes)
!-------------------------------------------------------------------------------
@ -195,13 +184,7 @@ module shapes
!-------------------------------------------------------------------------------
!
if (verbose) then
if (enabled) then
call print_parameter(verbose, "embedded shapes", "on" )
else
call print_parameter(verbose, "embedded shapes", "off")
end if
end if
if (verbose) call print_parameter(verbose, "embedded shapes", enabled, "on")
!-------------------------------------------------------------------------------
!

35
sources/sources.F90
View File

@ -31,6 +31,23 @@ module sources
implicit none
! interfaces for the extra source terms subroutine
!
abstract interface
subroutine update_extra_sources_iface(pdata, t, dt, du)
use blocks, only : block_data
implicit none
type(block_data), pointer , intent(inout) :: pdata
real(kind=8) , intent(in) :: t, dt
real(kind=8), dimension(:,:,:,:), intent(inout) :: du
end subroutine
end interface
! pointer to the extra source terms subroutine
!
procedure(update_extra_sources_iface), pointer, save :: &
update_extra_sources => null()
! GLM-MHD source terms type (1 - EGLM, 2 - HEGLM)
!
integer , save :: glm_type = 0
@ -46,14 +63,10 @@ module sources
real(kind=8) , save :: anomalous = 0.0d+00
real(kind=8) , save :: jcrit = 1.0d+00
! by default everything is private
!
private
! declare public subroutines
!
public :: initialize_sources, finalize_sources, print_sources
public :: update_sources
public :: update_sources, update_extra_sources
public :: viscosity, resistivity
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
@ -253,10 +266,10 @@ module sources
subroutine update_sources(pdata, t, dt, du)
use blocks , only : block_data
use coordinates , only : nn => bcells
use coordinates , only : ax, adx, ay, ady, adz
use coordinates, only : nn => bcells
use coordinates, only : ax, adx, ay, ady, adz
#if NDIMS == 3
use coordinates , only : az
use coordinates, only : az
#endif /* NDIMS == 3 */
use equations , only : magnetized
use equations , only : inx, iny, inz
@ -265,7 +278,6 @@ module sources
use gravity , only : gravity_enabled, gravitational_acceleration
use helpers , only : print_message
use operators , only : divergence, gradient, laplace, curl
use user_problem, only : update_user_sources
use workspace , only : resize_workspace, work, work_in_use
implicit none
@ -712,9 +724,10 @@ module sources
end if ! magnetized
! add user defined source terms
! add extra source terms
!
call update_user_sources(pdata, t, dt, du(:,:,:,:))
if (associated(update_extra_sources)) &
call update_extra_sources(pdata, t, dt, du(:,:,:,:))
100 continue

24
sources/system.F90
View File

@ -117,6 +117,7 @@ module system
use problems , only : initialize_problems
use random , only : reset_seeds
use statistics , only : initialize_statistics
use user_problem , only : initialize_user_problem
implicit none
@ -169,6 +170,11 @@ module system
call print_message(loc, "Could not initialize module PROBLEMS!")
if (check_status(status /= 0)) return
call initialize_user_problem(name, nrun, verbose, status)
if (status /=0) &
call print_message(loc, "Could not initialize module USER_PROBLEM!")
if (check_status(status /= 0)) return
call initialize_boundaries(status)
if (status /=0) &
call print_message(loc, "Could not initialize module BOUNDARIES!")
@ -217,6 +223,8 @@ module system
use mesh , only : finalize_mesh
use problems , only : finalize_problems
use statistics , only : finalize_statistics
use user_problem , only : finalize_user_problem
implicit none
@ -241,6 +249,10 @@ module system
if (status /=0) &
call print_message(loc, "Could not finalize module BOUNDARIES!")
call finalize_user_problem(status)
if (status /=0) &
call print_message(loc, "Could not finalize module USER_PROBLEM!")
call finalize_problems(status)
if (status /=0) &
call print_message(loc, "Could not finalize module PROBLEMS!")
@ -362,7 +374,11 @@ module system
return
end if
call boundary_variables(time, 0.0d+00)
call boundary_variables(time, 0.0d+00, status)
if (status /= 0) then
call print_message(loc, "Could not update variable boundaries!")
return
end if
else
@ -372,7 +388,11 @@ module system
return
end if
call boundary_variables(0.0d+00, 0.0d+00)
call boundary_variables(0.0d+00, 0.0d+00, status)
if (status /= 0) then
call print_message(loc, "Could not update variable boundaries!")
return
end if
call initialize_time_step()

13
sources/user_problem.F90
View File

@ -76,7 +76,7 @@ module user_problem
private
public :: initialize_user_problem, finalize_user_problem
public :: setup_user_problem, update_user_sources, update_user_shapes
public :: setup_user_problem
public :: user_boundary_x, user_boundary_y, user_boundary_z
public :: user_statistics
@ -103,6 +103,7 @@ module user_problem
!
subroutine initialize_user_problem(problem, rcount, verbose, status)
use boundaries , only : custom_boundary_x, custom_boundary_y
#if NDIMS == 3
use constants , only : pi
#endif /* NDIMS == 3 */
@ -110,8 +111,10 @@ module user_problem
use coordinates, only : ng => nghosts, ady, xlen, zlen, ymax
use equations , only : magnetized, csnd, csnd2
use helpers , only : print_section, print_parameter, print_message
use mesh , only : setup_problem
use parameters , only : get_parameter
use random , only : randuni, randsym
use shapes , only : update_shapes
implicit none
@ -404,6 +407,14 @@ module user_problem
call print_parameter(verbose, '( ) yabs' , yabs)
end if
! set procedure pointers for the problem setup subroutine, the shapes update
! for the absorption layer, and boundary conditions;
!
setup_problem => setup_user_problem
update_shapes => update_user_shapes
custom_boundary_x => user_boundary_x
custom_boundary_y => user_boundary_y
!-------------------------------------------------------------------------------
!
end subroutine initialize_user_problem