From 25ad944f444613f25ccc7573e7f7da4f5d3557be Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Fri, 26 Nov 2021 13:37:25 -0300 Subject: [PATCH 1/4] SOURCES: Add procedure interface for extra source terms. A user defined source terms can be pointer to the procedure pointer 'update_extra_sources'. Signed-off-by: Grzegorz Kowal --- sources/sources.F90 | 53 +++++++++++++++++++++++++--------------- sources/user_problem.F90 | 2 +- 2 files changed, 34 insertions(+), 21 deletions(-) diff --git a/sources/sources.F90 b/sources/sources.F90 index f0d1870..ada34c0 100644 --- a/sources/sources.F90 +++ b/sources/sources.F90 @@ -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 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - @@ -252,21 +265,20 @@ 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 blocks , only : block_data + 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 - use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ien - use equations , only : ibx, iby, ibz, ibp - 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 + use equations , only : magnetized + use equations , only : inx, iny, inz + use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ien + use equations , only : ibx, iby, ibz, ibp + use gravity , only : gravity_enabled, gravitational_acceleration + use helpers , only : print_message + use operators , only : divergence, gradient, laplace, curl + 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 diff --git a/sources/user_problem.F90 b/sources/user_problem.F90 index 3a9e75d..ca84d57 100644 --- a/sources/user_problem.F90 +++ b/sources/user_problem.F90 @@ -40,7 +40,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, update_user_shapes public :: user_boundary_x, user_boundary_y, user_boundary_z public :: user_statistics From 4ceb0089870f3ddc5356739ddb6fb72a4ed96bcb Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Fri, 26 Nov 2021 13:47:05 -0300 Subject: [PATCH 2/4] SHAPES: Do not pull update_user_shapes() from module USER_PROBLEM. The pointer to update_shapes can be associated in initialize_user_problem(). Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 2 +- sources/shapes.F90 | 25 ++++--------------------- sources/user_problem.F90 | 2 +- 3 files changed, 6 insertions(+), 23 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index b301be1..0a11315 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -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) !------------------------------------------------------------------------------- diff --git a/sources/shapes.F90 b/sources/shapes.F90 index 83f7a32..006c340 100644 --- a/sources/shapes.F90 +++ b/sources/shapes.F90 @@ -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") !------------------------------------------------------------------------------- ! diff --git a/sources/user_problem.F90 b/sources/user_problem.F90 index ca84d57..a271e32 100644 --- a/sources/user_problem.F90 +++ b/sources/user_problem.F90 @@ -40,7 +40,7 @@ module user_problem private public :: initialize_user_problem, finalize_user_problem - public :: setup_user_problem, update_user_shapes + public :: setup_user_problem public :: user_boundary_x, user_boundary_y, user_boundary_z public :: user_statistics From fd40fffc64ff13d5518ce8d0a37c484575260769 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Fri, 26 Nov 2021 16:42:53 -0300 Subject: [PATCH 3/4] BOUNDARIES: Add interface to custom boundaries. Also return the call status from block_boundary_specific(). Signed-off-by: Grzegorz Kowal --- sources/boundaries.F90 | 173 +++++++++++++++++++++-------------------- sources/evolution.F90 | 3 +- sources/system.F90 | 12 ++- 3 files changed, 101 insertions(+), 87 deletions(-) diff --git a/sources/boundaries.F90 b/sources/boundaries.F90 index 28c53a6..8a47d70 100644 --- a/sources/boundaries.F90 +++ b/sources/boundaries.F90 @@ -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 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - @@ -356,29 +372,23 @@ module boundaries ! ! Arguments: ! -! t, dt - time and time increment; +! 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 + 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 ! @@ -1369,14 +1379,13 @@ module boundaries ! ! Arguments: ! -! t, dt - time and time increment; +! 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 + 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 !------------------------------------------------------------------------------- ! @@ -5076,15 +5060,16 @@ module boundaries ! ! Arguments: ! -! nc - the direction; -! side - the side of the boundary; -! t, dt - time and time increment; -! x, y, z - the block coordinates; -! qn - the variable array; +! nc - the direction; +! side - the side of the boundary; +! 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 diff --git a/sources/evolution.F90 b/sources/evolution.F90 index 0a11315..3c50ef9 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -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)) diff --git a/sources/system.F90 b/sources/system.F90 index faf3440..5ddd4f5 100644 --- a/sources/system.F90 +++ b/sources/system.F90 @@ -362,7 +362,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 +376,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() From 8fb3abacbf4415ae493284333cdd0ee9e9d4a35e Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Fri, 26 Nov 2021 16:48:54 -0300 Subject: [PATCH 4/4] PROBLEM, SYSTEM: Move USER_PROBLEM initialization to SYSTEM. Signed-off-by: Grzegorz Kowal --- sources/problems.F90 | 8 -------- sources/system.F90 | 12 ++++++++++++ sources/user_problem.F90 | 1 + 3 files changed, 13 insertions(+), 8 deletions(-) diff --git a/sources/problems.F90 b/sources/problems.F90 index 9e2d851..8f82fe2 100644 --- a/sources/problems.F90 +++ b/sources/problems.F90 @@ -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 diff --git a/sources/system.F90 b/sources/system.F90 index 5ddd4f5..23b863f 100644 --- a/sources/system.F90 +++ b/sources/system.F90 @@ -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!") diff --git a/sources/user_problem.F90 b/sources/user_problem.F90 index a271e32..b4b0345 100644 --- a/sources/user_problem.F90 +++ b/sources/user_problem.F90 @@ -69,6 +69,7 @@ module user_problem use equations , only : eos, csnd, csnd2, adiabatic_index use helpers , only : print_section, print_parameter + use mesh , only : setup_problem use parameters, only : get_parameter implicit none