From 0bf427d5cc65927cc3e01629ee2992c707c250ad Mon Sep 17 00:00:00 2001 From: emorway-usgs Date: Wed, 21 Aug 2024 14:26:56 -0700 Subject: [PATCH 1/7] fix(GweInputData): large GWE models were throwing stackoverflow error --- src/Model/Connection/GweInterfaceModel.f90 | 3 - src/Model/GroundWaterEnergy/gwe.f90 | 4 - src/Model/ModelUtilities/GweInputData.f90 | 181 +++++++-------------- 3 files changed, 61 insertions(+), 127 deletions(-) diff --git a/src/Model/Connection/GweInterfaceModel.f90 b/src/Model/Connection/GweInterfaceModel.f90 index 1a6636dd362..ce22883260b 100644 --- a/src/Model/Connection/GweInterfaceModel.f90 +++ b/src/Model/Connection/GweInterfaceModel.f90 @@ -195,9 +195,6 @@ subroutine gweifmod_df(this) this%ia => this%dis%con%ia this%ja => this%dis%con%ja ! - ! -- Define shared data (cpw, rhow, latent heat of vaporization) - call this%gwecommon%gweshared_dat_df(this%neq) - ! ! -- Allocate model arrays, now that neq and nja are assigned call this%allocate_arrays() ! diff --git a/src/Model/GroundWaterEnergy/gwe.f90 b/src/Model/GroundWaterEnergy/gwe.f90 index 2d20666e7d3..947584b24c8 100644 --- a/src/Model/GroundWaterEnergy/gwe.f90 +++ b/src/Model/GroundWaterEnergy/gwe.f90 @@ -148,7 +148,6 @@ end subroutine gwe_cr subroutine gwe_df(this) ! -- modules use SimModule, only: store_error - use GweInputDataModule, only: gweshared_dat_df ! -- dummy class(GweModelType) :: this ! -- local @@ -181,9 +180,6 @@ subroutine gwe_df(this) this%ia => this%dis%con%ia this%ja => this%dis%con%ja ! - ! -- Define shared data (cpw, rhow, latent heat of vaporization) - call this%gwecommon%gweshared_dat_df(this%neq) - ! ! -- Allocate model arrays, now that neq and nja are assigned call this%allocate_arrays() ! diff --git a/src/Model/ModelUtilities/GweInputData.f90 b/src/Model/ModelUtilities/GweInputData.f90 index d5d302a4a5f..9914b4f917a 100644 --- a/src/Model/ModelUtilities/GweInputData.f90 +++ b/src/Model/ModelUtilities/GweInputData.f90 @@ -7,7 +7,6 @@ module GweInputDataModule private public :: GweInputDataType public :: gweshared_dat_cr - public :: gweshared_dat_df public :: set_gwe_dat_ptrs !> Data for sharing among multiple packages. Originally read in from @@ -30,78 +29,28 @@ module GweInputDataModule contains - ! -- public - procedure, public :: gweshared_dat_df + ! public procedure, public :: set_gwe_dat_ptrs procedure, public :: gweshared_dat_da - ! -- private - procedure, private :: allocate_shared_vars - procedure, private :: set_gwe_shared_scalars - procedure, private :: set_gwe_shared_arrays + ! private + procedure, private :: set_gwe_scalar_ptrs + procedure, private :: set_gwe_array_ptrs end type GweInputDataType contains -!> @brief Allocate the shared data -!< + !> @brief Allocate the shared data + !< subroutine gweshared_dat_cr(gwe_dat) - ! -- modules - ! -- dummy + ! modules + ! dummy type(GweInputDataType), pointer :: gwe_dat !< the input data block - ! - ! -- Create the object + + ! create the object allocate (gwe_dat) - ! - ! -- Return - return - end subroutine gweshared_dat_cr -!> @brief Define the shared data -!< - subroutine gweshared_dat_df(this, nodes) - ! -- dummy - class(GweInputDataType) :: this !< the input data block - integer(I4B), intent(in) :: nodes - ! -- local - ! - ! -- Allocate variables - call this%allocate_shared_vars(nodes) - ! - ! -- Return - return - end subroutine gweshared_dat_df - - !> @brief Define the information this object holds - !! - !! Allocate strings for storing label names - !! Intended to be analogous to allocate_scalars() - !< - subroutine allocate_shared_vars(this, nodes) - ! -- dummy - class(GweInputDataType) :: this !< TspLabelsType object - integer(I4B), intent(in) :: nodes - ! -- local - integer(I4B) :: i - ! - allocate (this%gwecpw) - allocate (this%gwerhow) - allocate (this%gwelatheatvap) - allocate (this%gwerhos(nodes)) - allocate (this%gwecps(nodes)) - ! - ! -- Initialize values - this%gwecpw = 4.184e3 !< J/(kg*C) @ 20C - this%gwerhow = DEP3 !< kg/m3 @ 20C - this%gwelatheatvap = 2.4535e6 !< J/kg @ 20C - do i = 1, nodes - this%gwecps(i) = DZERO - this%gwerhos(i) = DZERO - end do - ! - ! -- Return - return - end subroutine allocate_shared_vars + end subroutine gweshared_dat_cr !> @brief Allocate and read data from EST !! @@ -110,24 +59,23 @@ end subroutine allocate_shared_vars !! the aquifer material (cps), and density of the aquifer material !! (rhow) is used among other packages and is therefore stored in a !! separate class + !< subroutine set_gwe_dat_ptrs(this, gwerhow, gwecpw, gwelatheatvap, & gwerhos, gwecps) - ! -- dummy + ! dummy class(GweInputDataType) :: this !< the input data block - real(DP), intent(in) :: gwerhow !< ptr to density of water specified in EST - real(DP), intent(in) :: gwecpw !< ptr to heat capacity of water specified in EST - real(DP), intent(in) :: gwelatheatvap !< ptr to latent heat of vaporization specified in EST + real(DP), intent(in), pointer :: gwerhow !< ptr to density of water specified in EST + real(DP), intent(in), pointer :: gwecpw !< ptr to heat capacity of water specified in EST + real(DP), intent(in), pointer :: gwelatheatvap !< ptr to latent heat of vaporization specified in EST real(DP), dimension(:), pointer, contiguous :: gwerhos !< ptr to sptially-variably density of aquifer material specified in EST real(DP), dimension(:), pointer, contiguous :: gwecps !< ptr to sptially-variably heat capacity of aquifer material specified in EST - ! - ! -- Allocate scalars - call this%set_gwe_shared_scalars(gwerhow, gwecpw, gwelatheatvap) - ! - ! -- Allocate arrays - call this%set_gwe_shared_arrays(gwerhos, gwecps) - ! - ! -- Return - return + + ! allocate scalars + call this%set_gwe_scalar_ptrs(gwerhow, gwecpw, gwelatheatvap) + + ! allocate arrays + call this%set_gwe_array_ptrs(gwerhos, gwecps) + end subroutine set_gwe_dat_ptrs !> @brief Set pointers to scalars read by the EST package @@ -138,68 +86,61 @@ end subroutine set_gwe_dat_ptrs !! simulating evaporation will need access to latent heat of !! of vaporization. !< - subroutine set_gwe_shared_scalars(this, gwerhow, gwecpw, gwelatheatvap) - ! -- dummy + subroutine set_gwe_scalar_ptrs(this, gwerhow, gwecpw, gwelatheatvap) + ! dummy class(GweInputDataType) :: this !< GweInputDataType object - real(DP), intent(in) :: gwerhow - real(DP), intent(in) :: gwecpw - real(DP), intent(in), optional :: gwelatheatvap - ! - ! -- Set the pointers - ! -- Fixed density of water to be used by GWE - this%gwerhow = gwerhow - ! -- Spatially constant heat capacity of water ! kluge note: "specific heat" (which is heat capacity per unit mass) is probably the more correct term - this%gwecpw = gwecpw - ! -- Latent heat of vaporization + real(DP), pointer, intent(in) :: gwerhow !< density of water + real(DP), pointer, intent(in) :: gwecpw !< mass-based heat capacity of water + real(DP), pointer, intent(in), optional :: gwelatheatvap !< latent heat of vaporization + + ! set the pointers + ! fixed density of water to be used by GWE (vs GWT-based density) + this%gwerhow => gwerhow + ! spatially constant heat capacity of water ! kluge note: "specific heat" (which is heat capacity per unit mass) is probably the more correct term + this%gwecpw => gwecpw + ! latent heat of vaporization if (present(gwelatheatvap)) then - this%gwelatheatvap = gwelatheatvap + this%gwelatheatvap => gwelatheatvap end if - ! - ! -- Return - return - end subroutine set_gwe_shared_scalars + + end subroutine set_gwe_scalar_ptrs !> @brief Set pointers to data arrays read by the EST package !! for use by other packages !! - !! Set pointers to GWE-related arrays for use - !! by multiple packages. + !! Set pointers to GWE-related arrays for use by multiple packages !< - subroutine set_gwe_shared_arrays(this, gwerhos, gwecps) - ! -- dummy + subroutine set_gwe_array_ptrs(this, gwerhos, gwecps) + ! dummy class(GweInputDataType) :: this !< GweInputDataType object real(DP), dimension(:), pointer, contiguous, intent(in) :: gwerhos real(DP), dimension(:), pointer, contiguous, intent(in) :: gwecps - ! - ! -- Set the pointers - ! -- Spatially-variable density of aquifer solid material - this%gwerhos = gwerhos - ! -- Spatially-variable heat capacity of aquifer solid material - this%gwecps = gwecps - ! - ! -- Return - return - end subroutine set_gwe_shared_arrays + + ! set the pointers + ! spatially-variable density of aquifer solid material + this%gwerhos => gwerhos + ! spatially-variable heat capacity of aquifer solid material + this%gwecps => gwecps + + end subroutine set_gwe_array_ptrs !> @ brief Deallocate memory !! - !! Deallocate GWE shared data array memory + !! Set pointers to null !< subroutine gweshared_dat_da(this) - ! -- dummy + ! dummy class(GweInputDataType) :: this !< the input data block - ! - ! -- Scalars - deallocate (this%gwelatheatvap) - deallocate (this%gwerhow) - deallocate (this%gwecpw) - ! - ! -- Arrays - deallocate (this%gwerhos) - deallocate (this%gwecps) - ! - ! -- Return - return + + ! scalars + this%gwelatheatvap => null() + this%gwerhow => null() + this%gwecpw => null() + + ! arrays + this%gwerhos => null() + this%gwecps => null() + end subroutine gweshared_dat_da end module GweInputDataModule From 1496b6761513c65d34cb4ee997761ddd55794408 Mon Sep 17 00:00:00 2001 From: emorway-usgs Date: Mon, 26 Aug 2024 14:20:14 -0700 Subject: [PATCH 2/7] trying this fix a bit differently --- src/Model/Connection/GweInterfaceModel.f90 | 3 + src/Model/GroundWaterEnergy/gwe-est.f90 | 2 +- src/Model/GroundWaterEnergy/gwe.f90 | 4 + src/Model/ModelUtilities/GweInputData.f90 | 117 +++++++++++++++------ 4 files changed, 91 insertions(+), 35 deletions(-) diff --git a/src/Model/Connection/GweInterfaceModel.f90 b/src/Model/Connection/GweInterfaceModel.f90 index ce22883260b..1a6636dd362 100644 --- a/src/Model/Connection/GweInterfaceModel.f90 +++ b/src/Model/Connection/GweInterfaceModel.f90 @@ -195,6 +195,9 @@ subroutine gweifmod_df(this) this%ia => this%dis%con%ia this%ja => this%dis%con%ja ! + ! -- Define shared data (cpw, rhow, latent heat of vaporization) + call this%gwecommon%gweshared_dat_df(this%neq) + ! ! -- Allocate model arrays, now that neq and nja are assigned call this%allocate_arrays() ! diff --git a/src/Model/GroundWaterEnergy/gwe-est.f90 b/src/Model/GroundWaterEnergy/gwe-est.f90 index 6f0f8a9493f..2384da32209 100644 --- a/src/Model/GroundWaterEnergy/gwe-est.f90 +++ b/src/Model/GroundWaterEnergy/gwe-est.f90 @@ -150,7 +150,7 @@ subroutine est_ar(this, dis, ibound) ! ! -- set data required by other packages call this%gwecommon%set_gwe_dat_ptrs(this%rhow, this%cpw, this%latheatvap, & - this%rhos, this%cps) + this%rhos, this%cps, dis%nodes) ! ! -- Return return diff --git a/src/Model/GroundWaterEnergy/gwe.f90 b/src/Model/GroundWaterEnergy/gwe.f90 index 947584b24c8..2d20666e7d3 100644 --- a/src/Model/GroundWaterEnergy/gwe.f90 +++ b/src/Model/GroundWaterEnergy/gwe.f90 @@ -148,6 +148,7 @@ end subroutine gwe_cr subroutine gwe_df(this) ! -- modules use SimModule, only: store_error + use GweInputDataModule, only: gweshared_dat_df ! -- dummy class(GweModelType) :: this ! -- local @@ -180,6 +181,9 @@ subroutine gwe_df(this) this%ia => this%dis%con%ia this%ja => this%dis%con%ja ! + ! -- Define shared data (cpw, rhow, latent heat of vaporization) + call this%gwecommon%gweshared_dat_df(this%neq) + ! ! -- Allocate model arrays, now that neq and nja are assigned call this%allocate_arrays() ! diff --git a/src/Model/ModelUtilities/GweInputData.f90 b/src/Model/ModelUtilities/GweInputData.f90 index 9914b4f917a..ca3b830b0e4 100644 --- a/src/Model/ModelUtilities/GweInputData.f90 +++ b/src/Model/ModelUtilities/GweInputData.f90 @@ -7,6 +7,7 @@ module GweInputDataModule private public :: GweInputDataType public :: gweshared_dat_cr + public :: gweshared_dat_df public :: set_gwe_dat_ptrs !> Data for sharing among multiple packages. Originally read in from @@ -30,11 +31,13 @@ module GweInputDataModule contains ! public + procedure, public :: gweshared_dat_df procedure, public :: set_gwe_dat_ptrs procedure, public :: gweshared_dat_da ! private - procedure, private :: set_gwe_scalar_ptrs - procedure, private :: set_gwe_array_ptrs + procedure, private :: allocate_shared_vars + procedure, private :: set_gwe_shared_scalars + procedure, private :: set_gwe_shared_arrays end type GweInputDataType @@ -43,15 +46,55 @@ module GweInputDataModule !> @brief Allocate the shared data !< subroutine gweshared_dat_cr(gwe_dat) - ! modules ! dummy type(GweInputDataType), pointer :: gwe_dat !< the input data block - ! create the object + ! Create the object allocate (gwe_dat) end subroutine gweshared_dat_cr + !> @brief Define the shared data + !< + subroutine gweshared_dat_df(this, nodes) + ! dummy + class(GweInputDataType) :: this !< the input data block + integer(I4B), intent(in) :: nodes + + ! allocate variables + call this%allocate_shared_vars(nodes) + + end subroutine gweshared_dat_df + + !> @brief Define the information this object holds + !! + !! Allocate strings for storing label names + !! Intended to be analogous to allocate_scalars() + !< + subroutine allocate_shared_vars(this, nodes) + ! dummy + class(GweInputDataType) :: this !< TspLabelsType object + integer(I4B), intent(in) :: nodes + ! local + integer(I4B) :: i + + allocate (this%gwecpw) + allocate (this%gwerhow) + allocate (this%gwelatheatvap) + allocate (this%gwerhos(nodes)) + allocate (this%gwecps(nodes)) + + ! initialize values + this%gwecpw = 4.184e3 !< J/(kg*C) @ 20C + this%gwerhow = DEP3 !< kg/m3 @ 20C + this%gwelatheatvap = 2.4535e6 !< J/kg @ 20C + do i = 1, nodes + this%gwecps(i) = DZERO + this%gwerhos(i) = DZERO + end do + + end subroutine allocate_shared_vars + !> @brief Allocate and read data from EST !! !! EST data, including heat capacity of water (cpw), density of water @@ -61,20 +104,21 @@ end subroutine gweshared_dat_cr !! separate class !< subroutine set_gwe_dat_ptrs(this, gwerhow, gwecpw, gwelatheatvap, & - gwerhos, gwecps) + gwerhos, gwecps, nodes) ! dummy class(GweInputDataType) :: this !< the input data block - real(DP), intent(in), pointer :: gwerhow !< ptr to density of water specified in EST - real(DP), intent(in), pointer :: gwecpw !< ptr to heat capacity of water specified in EST - real(DP), intent(in), pointer :: gwelatheatvap !< ptr to latent heat of vaporization specified in EST + real(DP), intent(in) :: gwerhow !< ptr to density of water specified in EST + real(DP), intent(in) :: gwecpw !< ptr to heat capacity of water specified in EST + real(DP), intent(in) :: gwelatheatvap !< ptr to latent heat of vaporization specified in EST real(DP), dimension(:), pointer, contiguous :: gwerhos !< ptr to sptially-variably density of aquifer material specified in EST real(DP), dimension(:), pointer, contiguous :: gwecps !< ptr to sptially-variably heat capacity of aquifer material specified in EST + integer(I4B), intent(in) :: nodes ! allocate scalars - call this%set_gwe_scalar_ptrs(gwerhow, gwecpw, gwelatheatvap) + call this%set_gwe_shared_scalars(gwerhow, gwecpw, gwelatheatvap) ! allocate arrays - call this%set_gwe_array_ptrs(gwerhos, gwecps) + call this%set_gwe_shared_arrays(gwerhos, gwecps, nodes) end subroutine set_gwe_dat_ptrs @@ -86,60 +130,65 @@ end subroutine set_gwe_dat_ptrs !! simulating evaporation will need access to latent heat of !! of vaporization. !< - subroutine set_gwe_scalar_ptrs(this, gwerhow, gwecpw, gwelatheatvap) + subroutine set_gwe_shared_scalars(this, gwerhow, gwecpw, gwelatheatvap) ! dummy class(GweInputDataType) :: this !< GweInputDataType object - real(DP), pointer, intent(in) :: gwerhow !< density of water - real(DP), pointer, intent(in) :: gwecpw !< mass-based heat capacity of water - real(DP), pointer, intent(in), optional :: gwelatheatvap !< latent heat of vaporization + real(DP), intent(in) :: gwerhow + real(DP), intent(in) :: gwecpw + real(DP), intent(in), optional :: gwelatheatvap - ! set the pointers - ! fixed density of water to be used by GWE (vs GWT-based density) - this%gwerhow => gwerhow + ! fixed density of water to be used by GWE + this%gwerhow = gwerhow ! spatially constant heat capacity of water ! kluge note: "specific heat" (which is heat capacity per unit mass) is probably the more correct term - this%gwecpw => gwecpw + this%gwecpw = gwecpw ! latent heat of vaporization if (present(gwelatheatvap)) then - this%gwelatheatvap => gwelatheatvap + this%gwelatheatvap = gwelatheatvap end if - end subroutine set_gwe_scalar_ptrs + end subroutine set_gwe_shared_scalars !> @brief Set pointers to data arrays read by the EST package !! for use by other packages !! - !! Set pointers to GWE-related arrays for use by multiple packages + !! Set pointers to GWE-related arrays for use + !! by multiple packages. !< - subroutine set_gwe_array_ptrs(this, gwerhos, gwecps) + subroutine set_gwe_shared_arrays(this, gwerhos, gwecps, nodes) ! dummy class(GweInputDataType) :: this !< GweInputDataType object real(DP), dimension(:), pointer, contiguous, intent(in) :: gwerhos real(DP), dimension(:), pointer, contiguous, intent(in) :: gwecps + integer(I4B), intent(in) :: nodes + + ! local + integer(I4B) :: i - ! set the pointers - ! spatially-variable density of aquifer solid material - this%gwerhos => gwerhos - ! spatially-variable heat capacity of aquifer solid material - this%gwecps => gwecps + do i = 1, nodes + ! spatially-variable density of aquifer solid material + this%gwerhos(i) = gwerhos(i) + ! spatially-variable heat capacity of aquifer solid material + this%gwecps(i) = gwecps(i) + end do - end subroutine set_gwe_array_ptrs + end subroutine set_gwe_shared_arrays !> @ brief Deallocate memory !! - !! Set pointers to null + !! Deallocate GWE shared data array memory !< subroutine gweshared_dat_da(this) ! dummy class(GweInputDataType) :: this !< the input data block ! scalars - this%gwelatheatvap => null() - this%gwerhow => null() - this%gwecpw => null() + deallocate (this%gwelatheatvap) + deallocate (this%gwerhow) + deallocate (this%gwecpw) ! arrays - this%gwerhos => null() - this%gwecps => null() + deallocate (this%gwerhos) + deallocate (this%gwecps) end subroutine gweshared_dat_da From 7692f5269b1ee79647124037793541e4be926110 Mon Sep 17 00:00:00 2001 From: emorway-usgs Date: Tue, 27 Aug 2024 10:06:29 -0700 Subject: [PATCH 3/7] Back to 16cac6d with one additional change GweGweConnection.f90; more to follow. --- src/Model/Connection/GweGweConnection.f90 | 4 +- src/Model/Connection/GweInterfaceModel.f90 | 3 - src/Model/GroundWaterEnergy/gwe-est.f90 | 2 +- src/Model/GroundWaterEnergy/gwe.f90 | 4 - src/Model/ModelUtilities/GweInputData.f90 | 117 ++++++--------------- 5 files changed, 37 insertions(+), 93 deletions(-) diff --git a/src/Model/Connection/GweGweConnection.f90 b/src/Model/Connection/GweGweConnection.f90 index 2f4464b0d27..9abda79bf79 100644 --- a/src/Model/Connection/GweGweConnection.f90 +++ b/src/Model/Connection/GweGweConnection.f90 @@ -299,8 +299,8 @@ subroutine gwegwecon_ar(this) ! -- Set a pointer in the interface model to the gwecommon data if (this%gweModel%inest > 0) then - this%gweInterfaceModel%gwecommon%gwecpw = this%gweModel%gwecommon%gwecpw - this%gweInterfaceModel%gwecommon%gwerhow = this%gweModel%gwecommon%gwerhow + this%gweInterfaceModel%gwecommon%gwecpw => this%gweModel%gwecommon%gwecpw + this%gweInterfaceModel%gwecommon%gwerhow => this%gweModel%gwecommon%gwerhow end if !-- Set the equation scaling factor in the interface model to that of diff --git a/src/Model/Connection/GweInterfaceModel.f90 b/src/Model/Connection/GweInterfaceModel.f90 index 1a6636dd362..ce22883260b 100644 --- a/src/Model/Connection/GweInterfaceModel.f90 +++ b/src/Model/Connection/GweInterfaceModel.f90 @@ -195,9 +195,6 @@ subroutine gweifmod_df(this) this%ia => this%dis%con%ia this%ja => this%dis%con%ja ! - ! -- Define shared data (cpw, rhow, latent heat of vaporization) - call this%gwecommon%gweshared_dat_df(this%neq) - ! ! -- Allocate model arrays, now that neq and nja are assigned call this%allocate_arrays() ! diff --git a/src/Model/GroundWaterEnergy/gwe-est.f90 b/src/Model/GroundWaterEnergy/gwe-est.f90 index 2384da32209..6f0f8a9493f 100644 --- a/src/Model/GroundWaterEnergy/gwe-est.f90 +++ b/src/Model/GroundWaterEnergy/gwe-est.f90 @@ -150,7 +150,7 @@ subroutine est_ar(this, dis, ibound) ! ! -- set data required by other packages call this%gwecommon%set_gwe_dat_ptrs(this%rhow, this%cpw, this%latheatvap, & - this%rhos, this%cps, dis%nodes) + this%rhos, this%cps) ! ! -- Return return diff --git a/src/Model/GroundWaterEnergy/gwe.f90 b/src/Model/GroundWaterEnergy/gwe.f90 index 2d20666e7d3..947584b24c8 100644 --- a/src/Model/GroundWaterEnergy/gwe.f90 +++ b/src/Model/GroundWaterEnergy/gwe.f90 @@ -148,7 +148,6 @@ end subroutine gwe_cr subroutine gwe_df(this) ! -- modules use SimModule, only: store_error - use GweInputDataModule, only: gweshared_dat_df ! -- dummy class(GweModelType) :: this ! -- local @@ -181,9 +180,6 @@ subroutine gwe_df(this) this%ia => this%dis%con%ia this%ja => this%dis%con%ja ! - ! -- Define shared data (cpw, rhow, latent heat of vaporization) - call this%gwecommon%gweshared_dat_df(this%neq) - ! ! -- Allocate model arrays, now that neq and nja are assigned call this%allocate_arrays() ! diff --git a/src/Model/ModelUtilities/GweInputData.f90 b/src/Model/ModelUtilities/GweInputData.f90 index ca3b830b0e4..9914b4f917a 100644 --- a/src/Model/ModelUtilities/GweInputData.f90 +++ b/src/Model/ModelUtilities/GweInputData.f90 @@ -7,7 +7,6 @@ module GweInputDataModule private public :: GweInputDataType public :: gweshared_dat_cr - public :: gweshared_dat_df public :: set_gwe_dat_ptrs !> Data for sharing among multiple packages. Originally read in from @@ -31,13 +30,11 @@ module GweInputDataModule contains ! public - procedure, public :: gweshared_dat_df procedure, public :: set_gwe_dat_ptrs procedure, public :: gweshared_dat_da ! private - procedure, private :: allocate_shared_vars - procedure, private :: set_gwe_shared_scalars - procedure, private :: set_gwe_shared_arrays + procedure, private :: set_gwe_scalar_ptrs + procedure, private :: set_gwe_array_ptrs end type GweInputDataType @@ -46,55 +43,15 @@ module GweInputDataModule !> @brief Allocate the shared data !< subroutine gweshared_dat_cr(gwe_dat) + ! modules ! dummy type(GweInputDataType), pointer :: gwe_dat !< the input data block - ! Create the object + ! create the object allocate (gwe_dat) end subroutine gweshared_dat_cr - !> @brief Define the shared data - !< - subroutine gweshared_dat_df(this, nodes) - ! dummy - class(GweInputDataType) :: this !< the input data block - integer(I4B), intent(in) :: nodes - - ! allocate variables - call this%allocate_shared_vars(nodes) - - end subroutine gweshared_dat_df - - !> @brief Define the information this object holds - !! - !! Allocate strings for storing label names - !! Intended to be analogous to allocate_scalars() - !< - subroutine allocate_shared_vars(this, nodes) - ! dummy - class(GweInputDataType) :: this !< TspLabelsType object - integer(I4B), intent(in) :: nodes - ! local - integer(I4B) :: i - - allocate (this%gwecpw) - allocate (this%gwerhow) - allocate (this%gwelatheatvap) - allocate (this%gwerhos(nodes)) - allocate (this%gwecps(nodes)) - - ! initialize values - this%gwecpw = 4.184e3 !< J/(kg*C) @ 20C - this%gwerhow = DEP3 !< kg/m3 @ 20C - this%gwelatheatvap = 2.4535e6 !< J/kg @ 20C - do i = 1, nodes - this%gwecps(i) = DZERO - this%gwerhos(i) = DZERO - end do - - end subroutine allocate_shared_vars - !> @brief Allocate and read data from EST !! !! EST data, including heat capacity of water (cpw), density of water @@ -104,21 +61,20 @@ end subroutine allocate_shared_vars !! separate class !< subroutine set_gwe_dat_ptrs(this, gwerhow, gwecpw, gwelatheatvap, & - gwerhos, gwecps, nodes) + gwerhos, gwecps) ! dummy class(GweInputDataType) :: this !< the input data block - real(DP), intent(in) :: gwerhow !< ptr to density of water specified in EST - real(DP), intent(in) :: gwecpw !< ptr to heat capacity of water specified in EST - real(DP), intent(in) :: gwelatheatvap !< ptr to latent heat of vaporization specified in EST + real(DP), intent(in), pointer :: gwerhow !< ptr to density of water specified in EST + real(DP), intent(in), pointer :: gwecpw !< ptr to heat capacity of water specified in EST + real(DP), intent(in), pointer :: gwelatheatvap !< ptr to latent heat of vaporization specified in EST real(DP), dimension(:), pointer, contiguous :: gwerhos !< ptr to sptially-variably density of aquifer material specified in EST real(DP), dimension(:), pointer, contiguous :: gwecps !< ptr to sptially-variably heat capacity of aquifer material specified in EST - integer(I4B), intent(in) :: nodes ! allocate scalars - call this%set_gwe_shared_scalars(gwerhow, gwecpw, gwelatheatvap) + call this%set_gwe_scalar_ptrs(gwerhow, gwecpw, gwelatheatvap) ! allocate arrays - call this%set_gwe_shared_arrays(gwerhos, gwecps, nodes) + call this%set_gwe_array_ptrs(gwerhos, gwecps) end subroutine set_gwe_dat_ptrs @@ -130,65 +86,60 @@ end subroutine set_gwe_dat_ptrs !! simulating evaporation will need access to latent heat of !! of vaporization. !< - subroutine set_gwe_shared_scalars(this, gwerhow, gwecpw, gwelatheatvap) + subroutine set_gwe_scalar_ptrs(this, gwerhow, gwecpw, gwelatheatvap) ! dummy class(GweInputDataType) :: this !< GweInputDataType object - real(DP), intent(in) :: gwerhow - real(DP), intent(in) :: gwecpw - real(DP), intent(in), optional :: gwelatheatvap + real(DP), pointer, intent(in) :: gwerhow !< density of water + real(DP), pointer, intent(in) :: gwecpw !< mass-based heat capacity of water + real(DP), pointer, intent(in), optional :: gwelatheatvap !< latent heat of vaporization - ! fixed density of water to be used by GWE - this%gwerhow = gwerhow + ! set the pointers + ! fixed density of water to be used by GWE (vs GWT-based density) + this%gwerhow => gwerhow ! spatially constant heat capacity of water ! kluge note: "specific heat" (which is heat capacity per unit mass) is probably the more correct term - this%gwecpw = gwecpw + this%gwecpw => gwecpw ! latent heat of vaporization if (present(gwelatheatvap)) then - this%gwelatheatvap = gwelatheatvap + this%gwelatheatvap => gwelatheatvap end if - end subroutine set_gwe_shared_scalars + end subroutine set_gwe_scalar_ptrs !> @brief Set pointers to data arrays read by the EST package !! for use by other packages !! - !! Set pointers to GWE-related arrays for use - !! by multiple packages. + !! Set pointers to GWE-related arrays for use by multiple packages !< - subroutine set_gwe_shared_arrays(this, gwerhos, gwecps, nodes) + subroutine set_gwe_array_ptrs(this, gwerhos, gwecps) ! dummy class(GweInputDataType) :: this !< GweInputDataType object real(DP), dimension(:), pointer, contiguous, intent(in) :: gwerhos real(DP), dimension(:), pointer, contiguous, intent(in) :: gwecps - integer(I4B), intent(in) :: nodes - - ! local - integer(I4B) :: i - do i = 1, nodes - ! spatially-variable density of aquifer solid material - this%gwerhos(i) = gwerhos(i) - ! spatially-variable heat capacity of aquifer solid material - this%gwecps(i) = gwecps(i) - end do + ! set the pointers + ! spatially-variable density of aquifer solid material + this%gwerhos => gwerhos + ! spatially-variable heat capacity of aquifer solid material + this%gwecps => gwecps - end subroutine set_gwe_shared_arrays + end subroutine set_gwe_array_ptrs !> @ brief Deallocate memory !! - !! Deallocate GWE shared data array memory + !! Set pointers to null !< subroutine gweshared_dat_da(this) ! dummy class(GweInputDataType) :: this !< the input data block ! scalars - deallocate (this%gwelatheatvap) - deallocate (this%gwerhow) - deallocate (this%gwecpw) + this%gwelatheatvap => null() + this%gwerhow => null() + this%gwecpw => null() ! arrays - deallocate (this%gwerhos) - deallocate (this%gwecps) + this%gwerhos => null() + this%gwecps => null() end subroutine gweshared_dat_da From 490acb4d4a590d7f6c917ddb9b19355b7899d1e7 Mon Sep 17 00:00:00 2001 From: emorway-usgs Date: Wed, 28 Aug 2024 07:58:58 -0700 Subject: [PATCH 4/7] working on a new test that isn't quite ready, but needs to be under source control --- autotest/test_gwe_split_analyt.py | 6 +- src/Model/Connection/GweGweConnection.f90 | 287 ++++++++++++++-------- 2 files changed, 188 insertions(+), 105 deletions(-) diff --git a/autotest/test_gwe_split_analyt.py b/autotest/test_gwe_split_analyt.py index 8c77b03ad2b..132eb2d1b16 100644 --- a/autotest/test_gwe_split_analyt.py +++ b/autotest/test_gwe_split_analyt.py @@ -154,7 +154,7 @@ def calc_ener_input(primer_val): return ener_add_rate -# Define function to solve analytical solution +# Instatiate model to compare against analytical solution def assemble_half_model(sim, gwfname, gwfpath, side="right"): # Create GWF model gwf = flopy.mf6.MFModel( @@ -409,12 +409,12 @@ def build_models(idx, test): ) sim.register_ims_package(imsgwf, [gwf1.name, gwf2.name]) - # Create gw3 model + # Create first gwe model gwe1 = get_gwe_model( idx, sim, "energy1", "energy1", ener_input, side="left" ) - # Create gwe model + # Create second gwe model gwe2 = get_gwe_model( idx, sim, "energy2", "energy2", ener_input, side="right" ) diff --git a/src/Model/Connection/GweGweConnection.f90 b/src/Model/Connection/GweGweConnection.f90 index 9abda79bf79..1edfc2aba92 100644 --- a/src/Model/Connection/GweGweConnection.f90 +++ b/src/Model/Connection/GweGweConnection.f90 @@ -71,6 +71,7 @@ module GweGweConnectionModule procedure, private :: allocate_arrays procedure, private :: cfg_dist_vars procedure, private :: setGridExtent + procedure, private :: validateGweExchange procedure, private :: setFlowToExchange end type GweGweConnectionType @@ -80,37 +81,37 @@ module GweGweConnectionModule !> @brief Basic construction of the connection !< subroutine gweGweConnection_ctor(this, model, gweEx) - ! -- modules + ! modules use InputOutputModule, only: openfile - ! -- dummy + ! dummy class(GweGweConnectionType) :: this !< the connection class(NumericalModelType), pointer :: model !< the model owning this connection, !! this must be a GweModelType class(DisConnExchangeType), pointer :: gweEx !< the GWE-GWE exchange the interface model is created for - ! -- local + ! local character(len=LINELENGTH) :: fname character(len=LENCOMPONENTNAME) :: name class(*), pointer :: objPtr logical(LGP) :: write_ifmodel_listfile = .false. - ! + objPtr => model this%gweModel => CastAsGweModel(objPtr) objPtr => gweEx this%gweExchange => CastAsGweExchange(objPtr) - ! + if (gweEx%v_model1%is_local .and. gweEx%v_model2%is_local) then this%owns_exchange = associated(model, gweEx%model1) else this%owns_exchange = .true. end if - ! + if (gweEx%v_model1 == model) then write (name, '(a,i0)') 'GWECON1_', gweEx%id else write (name, '(a,i0)') 'GWECON2_', gweEx%id end if - ! - ! -- .lst file for interface model + + ! .lst file for interface model if (write_ifmodel_listfile) then fname = trim(name)//'.im.lst' call openfile(this%iout, 0, fname, 'LIST', filstat_opt='REPLACE') @@ -118,62 +119,62 @@ subroutine gweGweConnection_ctor(this, model, gweEx) trim(this%gweModel%name), 'from exchange ', & trim(gweEx%name) end if - ! - ! -- First call base constructor + + ! first call base constructor call this%SpatialModelConnectionType%spatialConnection_ctor(model, & gweEx, & name) - ! + call this%allocate_scalars() this%typename = 'GWE-GWE' this%iIfaceAdvScheme = 0 this%iIfaceXt3d = 0 this%exgflowSign = 1 - ! + allocate (this%gweInterfaceModel) this%interface_model => this%gweInterfaceModel - ! + end subroutine gweGweConnection_ctor !> @brief Allocate scalar variables for this connection !< subroutine allocate_scalars(this) - ! -- dummy + ! dummy class(GweGweConnectionType) :: this !< the connection - ! + call mem_allocate(this%iIfaceAdvScheme, 'IADVSCHEME', this%memoryPath) call mem_allocate(this%iIfaceXt3d, 'IXT3D', this%memoryPath) call mem_allocate(this%exgflowSign, 'EXGFLOWSIGN', this%memoryPath) - ! + end subroutine allocate_scalars !> @brief define the GWE-GWE connection !< subroutine gwegwecon_df(this) - ! -- dummy + ! dummy class(GweGweConnectionType) :: this !< the connection - ! -- local + ! local character(len=LENCOMPONENTNAME) :: imName - ! -- Determine advection scheme (the GWE-GWE exchange + ! determine advection scheme (the GWE-GWE exchange ! has been read at this point) this%iIfaceAdvScheme = this%gweExchange%iAdvScheme - ! - ! -- Determine xt3d setting on interface + + ! determine xt3d setting on interface this%iIfaceXt3d = this%gweExchange%ixt3d - ! -- Turn off when off in the owning model + ! turn off when off in the owning model if (this%gweModel%incnd > 0) then this%iIfaceXt3d = this%gweModel%cnd%ixt3d end if - ! -- Determine the required size of the interface model grid + ! determine the required size of the interface model grid call this%setGridExtent() - ! -- Now set up the GridConnection + ! now set up the GridConnection call this%spatialcon_df() - ! -- We have to 'catch up' and create the interface model + ! we have to 'catch up' and create the interface model ! here, then the remainder of this routine will be define if (this%prim_exchange%v_model1 == this%owner) then write (imName, '(a,i0)') 'GWEIM1_', this%gweExchange%id @@ -193,14 +194,14 @@ subroutine gwegwecon_df(this) call this%allocate_arrays() call this%gweInterfaceModel%allocate_fmi() - ! -- Connect X, RHS, IBOUND, and flowja + ! connect X, RHS, IBOUND, and flowja call this%spatialcon_setmodelptrs() - ! -- Connect pointers (used by BUY) + ! connect pointers (used by BUY) this%conc => this%gweInterfaceModel%x this%icbound => this%gweInterfaceModel%ibound - ! -- Add connections from the interface model to solution matrix + ! add connections from the interface model to solution matrix call this%spatialcon_connect() end subroutine gwegwecon_df @@ -208,16 +209,16 @@ end subroutine gwegwecon_df !> @brief Configure distributed variables for this interface model !< subroutine cfg_dist_vars(this) - ! -- dummy + ! dummy class(GweGweConnectionType) :: this !< the connection - ! + call this%cfg_dv('X', '', SYNC_NDS, & (/STG_BFR_CON_AR, STG_BFR_EXG_AD, STG_BFR_EXG_CF/)) call this%cfg_dv('IBOUND', '', SYNC_NDS, (/STG_BFR_CON_AR/)) call this%cfg_dv('TOP', 'DIS', SYNC_NDS, (/STG_BFR_CON_AR/)) call this%cfg_dv('BOT', 'DIS', SYNC_NDS, (/STG_BFR_CON_AR/)) call this%cfg_dv('AREA', 'DIS', SYNC_NDS, (/STG_BFR_CON_AR/)) - ! + if (this%gweInterfaceModel%cnd%idisp > 0) then call this%cfg_dv('ALH', 'CND', SYNC_NDS, (/STG_BFR_CON_AR/)) call this%cfg_dv('ALV', 'CND', SYNC_NDS, (/STG_BFR_CON_AR/)) @@ -233,11 +234,11 @@ subroutine cfg_dist_vars(this) call this%cfg_dv('GWFFLOWJA', 'FMI', SYNC_CON, (/STG_BFR_EXG_AD/)) call this%cfg_dv('GWFFLOWJA', 'FMI', SYNC_EXG, (/STG_BFR_EXG_AD/), & exg_var_name='GWFSIMVALS') - ! -- Fill porosity from est packages, needed for cnd + ! fill porosity from est packages, needed for cnd if (this%gweModel%incnd > 0 .and. this%gweModel%inest > 0) then call this%cfg_dv('POROSITY', 'EST', SYNC_NDS, (/STG_AFT_CON_AR/)) end if - ! + end subroutine cfg_dist_vars !> @brief Allocate array variables for this connection @@ -253,14 +254,14 @@ end subroutine allocate_arrays !> @brief Set required extent of the interface grid from !< the configuration subroutine setGridExtent(this) - ! -- dummy + ! dummy class(GweGweConnectionType) :: this !< the connection - ! -- local + ! local logical(LGP) :: hasAdv, hasCnd - ! + hasAdv = this%gweModel%inadv > 0 hasCnd = this%gweModel%incnd > 0 - ! + if (hasAdv) then if (this%iIfaceAdvScheme == 2) then this%exg_stencil_depth = 2 @@ -269,7 +270,7 @@ subroutine setGridExtent(this) end if end if end if - ! + if (hasCnd) then if (this%iIfaceXt3d > 0) then this%exg_stencil_depth = 2 @@ -278,7 +279,7 @@ subroutine setGridExtent(this) end if end if end if - ! + end subroutine setGridExtent !> @brief allocate and read/set the connection's data structures @@ -291,25 +292,25 @@ subroutine gwegwecon_ar(this) ! called, which is why we do it here call this%validateConnection() - ! -- Allocate and read base + ! allocate and read base call this%spatialcon_ar() ! ... and now the interface model call this%gweInterfaceModel%model_ar() - ! -- Set a pointer in the interface model to the gwecommon data + ! set a pointer in the interface model to the gwecommon data if (this%gweModel%inest > 0) then this%gweInterfaceModel%gwecommon%gwecpw => this%gweModel%gwecommon%gwecpw this%gweInterfaceModel%gwecommon%gwerhow => this%gweModel%gwecommon%gwerhow end if - !-- Set the equation scaling factor in the interface model to that of + ! set the equation scaling factor in the interface model to that of ! underlying GWE model if (this%gweModel%incnd > 0) then this%gweInterfaceModel%ieqnsclfac = this%gweModel%cnd%eqnsclfac end if - ! -- AR the movers and obs through the exchange + ! AR the movers and obs through the exchange if (this%owns_exchange) then !cdl implement this when MVT is ready !cdl if (this%gweExchange%inmvt > 0) then @@ -329,14 +330,17 @@ subroutine validateConnection(this) use SimModule, only: count_errors, store_error class(GweGweConnectionType) :: this !< this connection - ! -- Base validation, the spatial/geometry part + ! base validation, the spatial/geometry part call this%SpatialModelConnectionType%validateConnection() - ! -- We cannot validate this (yet) in parallel mode + ! we cannot validate this (yet) in parallel mode if (.not. this%gweExchange%v_model1%is_local) return if (.not. this%gweExchange%v_model2%is_local) return + + ! check specific cross-interface options/values that should be the same + call this%validateGweExchange() - ! -- GWE related matters + ! GWE related matters if ((this%gweExchange%gwemodel1%inadv > 0 .and. & this%gweExchange%gwemodel2%inadv == 0) .or. & (this%gweExchange%gwemodel2%inadv > 0 .and. & @@ -346,7 +350,7 @@ subroutine validateConnection(this) &and the other one is not' call store_error(errmsg) end if - ! + if ((this%gweExchange%gwemodel1%incnd > 0 .and. & this%gweExchange%gwemodel2%incnd == 0) .or. & (this%gweExchange%gwemodel2%incnd > 0 .and. & @@ -356,180 +360,259 @@ subroutine validateConnection(this) &and the other one is not' call store_error(errmsg) end if - ! - ! Abort on errors + + ! abort on errors if (count_errors() > 0) then write (errmsg, '(a)') 'Errors occurred while processing exchange(s)' call ustop() end if - ! + end subroutine validateConnection subroutine gwegwecon_rp(this) - ! -- dummy + ! dummy class(GweGweConnectionType) :: this !< the connection - ! - ! Call exchange rp routines + + ! call exchange rp routines if (this%owns_exchange) then call this%gweExchange%exg_rp() end if - ! + end subroutine gwegwecon_rp !> @brief Advance this connection !< subroutine gwegwecon_ad(this) - ! + class(GweGweConnectionType) :: this !< this connection - ! - ! -- Recalculate conduction ellipse + + ! recalculate conduction ellipse if (this%gweInterfaceModel%incnd > 0) call this%gweInterfaceModel%cnd%cnd_ad() - ! + if (this%owns_exchange) then call this%gweExchange%exg_ad() end if - ! + end subroutine gwegwecon_ad subroutine gwegwecon_fc(this, kiter, matrix_sln, rhs_sln, inwtflag) - ! -- dummy + ! dummy class(GweGweConnectionType) :: this !< the connection integer(I4B), intent(in) :: kiter !< the iteration counter class(MatrixBaseType), pointer :: matrix_sln !< the system matrix real(DP), dimension(:), intent(inout) :: rhs_sln !< global right-hand-side integer(I4B), optional, intent(in) :: inwtflag !< newton-raphson flag - ! + call this%SpatialModelConnectionType%spatialcon_fc( & kiter, matrix_sln, rhs_sln, inwtflag) - ! - ! FC the movers through the exchange + + ! _fc the movers through the exchange if (this%owns_exchange) then if (this%gweExchange%inmvt > 0) then call this%gweExchange%mvt%mvt_fc(this%gweExchange%gwemodel1%x, & this%gweExchange%gwemodel2%x) end if end if - ! + end subroutine gwegwecon_fc subroutine gwegwecon_cq(this, icnvg, isuppress_output, isolnid) - ! -- dummy + ! dummy class(GweGweConnectionType) :: this !< the connection integer(I4B), intent(inout) :: icnvg !< convergence flag integer(I4B), intent(in) :: isuppress_output !< suppress output when =1 integer(I4B), intent(in) :: isolnid !< solution id - ! + call this%gweInterfaceModel%model_cq(icnvg, isuppress_output) call this%setFlowToExchange() - ! + end subroutine gwegwecon_cq !> @brief Set the flows (flowja from interface model) to the !< simvals in the exchange, leaving the budget calcution in there subroutine setFlowToExchange(this) - ! -- modules + ! modules use IndexMapModule - ! -- dummy + ! dummy class(GweGweConnectionType) :: this !< this connection - ! -- local + ! local integer(I4B) :: i class(GweExchangeType), pointer :: gweEx type(IndexMapSgnType), pointer :: map - ! + if (this%owns_exchange) then gweEx => this%gweExchange map => this%interface_map%exchange_maps(this%interface_map%prim_exg_idx) - ! - ! -- Use (half of) the exchange map in reverse: + + ! use (half of) the exchange map in reverse: do i = 1, size(map%src_idx) if (map%sign(i) < 0) cycle ! simvals is defined from exg%m1 => exg%m2 gweEx%simvals(map%src_idx(i)) = & this%gweInterfaceModel%flowja(map%tgt_idx(i)) end do end if - ! + end subroutine setFlowToExchange subroutine gwegwecon_bd(this, icnvg, isuppress_output, isolnid) - ! -- modules + ! modules use BudgetModule, only: rate_accumulator - ! -- dummy + ! dummy class(GweGweConnectionType) :: this !< the connection integer(I4B), intent(inout) :: icnvg !< convergence flag integer(I4B), intent(in) :: isuppress_output !< suppress output when =1 integer(I4B), intent(in) :: isolnid !< solution id - ! - ! -- Call exchange budget routine, also calls bd + + ! call exchange budget routine, also calls _bd ! for movers. if (this%owns_exchange) then call this%gweExchange%exg_bd(icnvg, isuppress_output, isolnid) end if - ! + end subroutine gwegwecon_bd subroutine gwegwecon_ot(this) - ! -- dummy + ! dummy class(GweGweConnectionType) :: this !< the connection - ! - ! -- Call exg_ot() here as it handles all output processing + + ! call exg_ot() here as it handles all output processing ! based on gweExchange%simvals(:), which was correctly ! filled from gwegwecon if (this%owns_exchange) then call this%gweExchange%exg_ot() end if - ! + end subroutine gwegwecon_ot + !> @brief Validate the exchange, intercepting those + !! cases where two models have to be connected with an interface + !! model, where the individual configurations don't allow this + !! + !! Stops with error message on config mismatch + !< + subroutine validateGweExchange(this) + ! modules + use SimVariablesModule, only: errmsg, simulation_mode + use SimModule, only: store_error + use GweEstModule, only: GweEstType + + ! dummy + class(GweGweConnectionType) :: this !< this connection + + ! local + class(GweExchangeType), pointer :: gweEx + class(*), pointer :: modelPtr + class(GweModelType), pointer :: gweModel1 + class(GweModelType), pointer :: gweModel2 + type(GweEstType), pointer :: est1, est2 + logical(LGP) :: compatible + + gweEx => this%gweExchange + + ! we cannot validate the remainder (yet) in parallel mode + if (.not. gweEx%v_model1%is_local) return + if (.not. gweEx%v_model2%is_local) return + + modelPtr => this%gweExchange%model1 + gweModel1 => CastAsGweModel(modelPtr) + modelPtr => this%gweExchange%model2 + gweModel2 => CastAsGweModel(modelPtr) + + ! check that EST package usage is the same on both side of the interface + if ((gweModel1%inest > 0 .and. gweModel2%inest == 0) .or. & + (gweModel1%inest == 0 .and. gweModel2%inest > 0)) then + write (errmsg, '(2a)') 'Energy Storage and Transfer package should '// & + 'be enabled/disabled simultaneously in models connected with the '// & + 'interface model for exchange ', & + trim(gweEx%name) + call store_error(errmsg) + end if + + ! conduction options need to be the same in both model + if ((gweModel1%cnd%ixt3d > 0 .and. gweModel2%cnd%ixt3d == 0) .or. & + (gweModel1%cnd%ixt3d == 0 .and. gweModel2%cnd%ixt3d > 0)) then + write (errmsg, '(2a)') 'Use of XT3D to calculate conduction should '// & + 'be the same in both models, either both use XT3D or neither for '// & + ' exchange '// trim(gweEx%name) + call store_error(errmsg) + end if + + ! check compatibility of Energy Storage and Transfer (EST) package + compatible = .true. + est1 => gweModel1%est + est2 => gweModel2%est + if (est1%rhow /= est2%rhow) compatible = .false. + if (est1%cpw /= est2%cpw) compatible = .false. + ! if (est1%nrhospecies /= est2%nrhospecies) compatible = .false. + ! if (.not. all(buy1%drhodc == buy2%drhodc)) compatible = .false. + ! if (.not. all(buy1%crhoref == buy2%crhoref)) compatible = .false. + !if (.not. all(buy1%cauxspeciesname == buy2%cauxspeciesname)) then + ! compatible = .false. + !end if + if (.not. compatible) then + write (errmsg, '(6a)') 'Energy storage and transfer (EST) packages ', & + 'in model '// trim(gweEx%model1%name), ' and ', & + trim(gweEx%model2%name), & + ' should be equivalent to construct an '// & + ' interface model for exchange ', & + trim(gweEx%name) + call store_error(errmsg) + end if + + end subroutine validateGweExchange + + !> @brief Deallocate all resources + !< subroutine gwegwecon_da(this) - ! -- dummy + ! dummy class(GweGweConnectionType) :: this !< the connection - ! -- local + ! local logical(LGP) :: isOpen - ! - ! -- Scalars + + ! scalars call mem_deallocate(this%iIfaceAdvScheme) call mem_deallocate(this%iIfaceXt3d) call mem_deallocate(this%exgflowSign) - ! - ! -- Arrays + + ! arrays call mem_deallocate(this%exgflowjaGwe) - ! - ! -- Interface model + + ! interface model call this%gweInterfaceModel%model_da() deallocate (this%gweInterfaceModel) - ! - ! -- Dealloc base + + ! dealloc base call this%spatialcon_da() - ! + inquire (this%iout, opened=isOpen) if (isOpen) then close (this%iout) end if - ! - ! -- We need to deallocate the exchange we own: + + ! we need to deallocate the exchange we own: if (this%owns_exchange) then call this%gweExchange%exg_da() end if - ! + end subroutine gwegwecon_da !> @brief Cast to GweGweConnectionType !< function CastAsGweGweConnection(obj) result(res) implicit none - ! -- dummy + ! dummy class(*), pointer, intent(inout) :: obj !< object to be cast - ! -- return + ! return class(GweGweConnectionType), pointer :: res !< the GweGweConnection - ! + res => null() if (.not. associated(obj)) return - ! + select type (obj) class is (GweGweConnectionType) res => obj end select - ! + end function CastAsGweGweConnection end module From d1ba03b6ce758f8e99bd51db62a16a6fa98cf788 Mon Sep 17 00:00:00 2001 From: emorway-usgs Date: Wed, 28 Aug 2024 14:34:28 -0700 Subject: [PATCH 5/7] add new test --- autotest/test_gwe_bad_input.py | 507 +++++++++++++++++++++++++++++++++ 1 file changed, 507 insertions(+) create mode 100644 autotest/test_gwe_bad_input.py diff --git a/autotest/test_gwe_bad_input.py b/autotest/test_gwe_bad_input.py new file mode 100644 index 00000000000..2db94c5ebe8 --- /dev/null +++ b/autotest/test_gwe_bad_input.py @@ -0,0 +1,507 @@ +""" +Test problem for GWE designed to fail when neighboring gwe models +have different values of DENSITY_WATER specified. For this case, it's +not clear what value should be used in the interface model. As such, +model should exit with error msg informing the user of the issue. + +Problem set up: + + "Left" Model "Right" Model + ------------ ------------- + + +-------+-------+-------+ +-------+-------+-------+ + | | | | | | | | + | | | | <- connected -> | | | | + | | | | | | | | + +-------+-------+-------+ +-------+-------+-------+ + | | | | | | | | + | | | | <- connected -> | | | | + | | | | | | | | + +-------+-------+-------+ +-------+-------+-------+ + + ---> Direction of flow ---> +""" + +# Imports + +import os + +import flopy +import numpy as np +import pytest +from framework import TestFramework + +# Base simulation and model name and workspace + +scheme = "UPSTREAM" +# scheme = "TVD" + +cases = ["good_spec", "bad_spec"] +rhow = [999.0, 1001.0] + +# Model units +length_units = "meters" +time_units = "days" + +# Table MODFLOW 6 GWF & GWE params + +nrow = 1 +ncol = 3 +nlay = 2 +delr = 10.0 # Column width ($m$) +delc = 10.0 # Row width ($m$) +top = 20.0 # Top of model +k11 = 100.0 # Horizontal hydraulic conductivity ($m/d$) +ss = 1e-6 # Specific storage +sy = 0.20 # Specific Yield +prsity = 0.20 # Porosity +nper = 2 # Number of periods +perlen = [1, 100] # Simulation time ($days$) +nstp = [1, 10] # 10 day transient time steps +steady = {0: True, 1: False} +transient = {0: False, 1: True} + +# Set some static model parameter values + +k33 = k11 # Vertical hydraulic conductivity ($m/d$) +idomain = 1 # All cells included in the simulation +iconvert = 1 # All cells are convertible + +icelltype = 1 # Cell conversion type (>1: unconfined) + +# Set some static transport related model parameter values +botm = [] +botm.append(np.ones((nrow, ncol), dtype=float) * 10) +botm.append(np.zeros((nrow, ncol), dtype=float)) +botm = np.array(botm) + +# GWE related parameters +cpw = 4183.0 +lhv = 2454.0 +cps = 760.0 +rhos = 1500.0 + +# Head input +left_hd = 15.0 +right_hd = 2.0 +strt_hd1 = np.ones((nrow, ncol), dtype=float) * 11.0 +strt_hd2 = np.ones((nrow, ncol), dtype=float) * 11.0 +strt_hd1[0] = strt_hd2[0] = left_hd +strt_hd1[-1] = strt_hd2[-1] = right_hd +strt_hd = np.array([strt_hd1, strt_hd2]) +strt_temp = 10.0 + +chd_left = {} +chd_dat_left = [ + [(0, 0, 0), left_hd], + [(1, 0, 0), left_hd], +] +chd_left.update({0: chd_dat_left}) + +chd_right = {} +(chd_right.update({0: [[(1, 0, ncol - 1), right_hd]]}),) + +ctp_left = {} +ctp_dat_left = [ + [(0, 0, 0), strt_temp], + [(1, 0, 0), strt_temp], +] +ctp_left.update({0: ctp_dat_left}) + + +dispersivity = 0.0 # dispersion (remember, 1D model) + +# Set solver parameter values (and related) +nouter, ninner = 100, 300 +hclose, rclose, relax = 1e-10, 1e-10, 1.0 +ttsmult = 1.0 + +# Set up temporal data used by TDIS file +tdis_rc = [] +for i in np.arange(nper): + tdis_rc.append((perlen[i], nstp[i], ttsmult)) + +# ### Generate MODFLOW 6 Example test model +# + + +def build_gwf_model(sim, gwfname, side="left"): + pckg_suffix = "-1" + if side == "right": + pckg_suffix = "-2" + + # Instantiating MODFLOW 6 groundwater flow model + gwf = flopy.mf6.ModflowGwf( + sim, + modelname=gwfname, + newtonoptions="UNDER_RELAXATION", + save_flows=True, + model_nam_file=f"{gwfname}.nam", + ) + gwf.name_file.save_flows = True + + xorigin = 0.0 + if side == "right": + xorigin = ncol * delr + + # Instantiating MODFLOW 6 discretization package + flopy.mf6.ModflowGwfdis( + gwf, + xorigin=xorigin, + length_units=length_units, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + idomain=np.ones((nlay, nrow, ncol)), + pname="DIS" + pckg_suffix, + filename=f"{gwfname}.dis", + ) + + # Instantiating MODFLOW 6 storage package + flopy.mf6.ModflowGwfsto( + gwf, + ss=ss, + sy=sy, + iconvert=iconvert, + steady_state=steady, + transient=transient, + pname="STO" + pckg_suffix, + filename=f"{gwfname}.sto", + ) + + # Instantiating MODFLOW 6 node-property flow package + flopy.mf6.ModflowGwfnpf( + gwf, + save_flows=True, + icelltype=icelltype, + k=k11, + k33=k33, + save_specific_discharge=True, + pname="NPF" + pckg_suffix, + filename=f"{gwfname}.npf", + ) + + # Instantiating MODFLOW 6 constant head package + if side == "left": + chd_mf6 = chd_left + elif side == "right": + chd_mf6 = chd_right + + flopy.mf6.ModflowGwfchd( + gwf, + stress_period_data=chd_mf6, + pname="CHD" + pckg_suffix, + filename=f"{gwfname}.chd", + ) + + # Instantiating MODFLOW 6 initial conditions package for flow model + flopy.mf6.ModflowGwfic( + gwf, strt=strt_hd, pname="IC" + pckg_suffix, filename=f"{gwfname}.ic" + ) + + # Instantiating MODFLOW 6 output control package for flow model + flopy.mf6.ModflowGwfoc( + gwf, + pname="OC" + pckg_suffix, + head_filerecord=f"{gwfname}.hds", + budget_filerecord=f"{gwfname}.cbc", + headprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")], + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + ) + + gwf.set_model_relative_path(gwfname) + return gwf + + +# Build GWE model +def build_gwe_model(idx, sim, gwename, side="left"): + # Set package suffix + pckg_suffix = "-1" + if side == "right": + pckg_suffix = "-2" + + # Set origin information + xorigin = 0.0 + if side == "right": + xorigin = ncol * delr + + # Instantiate GWE model + gwe = flopy.mf6.ModflowGwe( + sim, modelname=gwename, model_nam_file=f"{gwename}.nam" + ) + gwe.name_file.save_flows = True + + # Instantiating MODFLOW 6 transport discretization package + flopy.mf6.ModflowGwedis( + gwe, + xorigin=xorigin, + length_units=length_units, + nogrb=True, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + idomain=1, + pname="DIS" + pckg_suffix, + filename=f"{gwename}.dis", + ) + + # Instantiating MODFLOW 6 transport initial concentrations + flopy.mf6.ModflowGweic( + gwe, strt=strt_temp, pname="IC" + pckg_suffix, filename=f"{gwename}.ic" + ) + + # Instantiating MODFLOW 6 transport advection package + flopy.mf6.ModflowGweadv( + gwe, + scheme=scheme, + pname="ADV" + pckg_suffix, + filename=f"{gwename}.adv", + ) + + # Instantiating MODFLOW 6 transport dispersion package + flopy.mf6.ModflowGwecnd( + gwe, + xt3d_off=False, + alh=dispersivity, + ath1=dispersivity, + ktw=0.5918 * 86400, + kts=0.2700 * 86400, + pname="CND" + pckg_suffix, + filename=f"{gwename}.cnd", + ) + + # Instantiating MODFLOW 6 transport mass storage package (formerly "reaction" package in MT3DMS) + rhow_mf6 = rhow[0] + if idx > 0 and side == "right": + # Set a parameter value that should trip a failure + rhow_mf6 = rhow[idx] + + flopy.mf6.ModflowGweest( + gwe, + save_flows=True, + porosity=prsity, + heat_capacity_water=cpw, + density_water=rhow_mf6, + latent_heat_vaporization=lhv, + heat_capacity_solid=cps, + density_solid=rhos, + pname="EST" + pckg_suffix, + filename=f"{gwename}.est", + ) + + # Instantiating MODFLOW 6 heat transport source-sink mixing package + flopy.mf6.ModflowGwessm( + gwe, sources=[[]], pname="SSM" + pckg_suffix, filename=f"{gwename}.ssm" + ) + + # Instantiate MODFLOW 6 heat transport output control package + flopy.mf6.ModflowGweoc( + gwe, + pname="OC" + pckg_suffix, + budget_filerecord=f"{gwename}.cbc", + temperature_filerecord=f"{gwename}.ucn", + temperatureprintrecord=[ + ("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL") + ], + saverecord=[("TEMPERATURE", "ALL"), ("BUDGET", "ALL")], + printrecord=[("TEMPERATURE", "ALL"), ("BUDGET", "ALL")], + ) + + # Instantiate a constant temperature in 1 of the dry cells + if side == "left": + flopy.mf6.ModflowGwectp( + gwe, + stress_period_data=ctp_left, + pname="CTP" + pckg_suffix, + filename=f"{gwename}.ctp", + ) + + return gwe + + +def build_models(idx, test): + # Base MF6 GWF model type + ws = test.workspace + name = cases[idx] + + print(f"Building MF6 model...{name}") + + # generate names for each model + gwename1 = "gwe-" + name + + sim = flopy.mf6.MFSimulation( + sim_name=name, sim_ws=ws, exe_name="mf6", version="mf6" + ) + + # Instantiating MODFLOW 6 time discretization + flopy.mf6.ModflowTdis( + sim, nper=nper, perioddata=tdis_rc, time_units=time_units + ) + + # left model + gwf1 = build_gwf_model(sim, "gwfleft", side="left") + + # right model + gwf2 = build_gwf_model(sim, "gwfright", side="right") + + # Add the exchange data + exgdata = [ + ((0, 0, ncol - 1), (0, 0, 0), 1, delr / 2, delr / 2, delc, 0.0, delr), + ((1, 0, ncol - 1), (1, 0, 0), 1, delr / 2, delr / 2, delc, 0.0, delr), + ] + flopy.mf6.ModflowGwfgwf( + sim, + exgtype="GWF6-GWF6", + nexg=len(exgdata), + exgmnamea=gwf1.name, + exgmnameb=gwf2.name, + exchangedata=exgdata, + xt3d=True, + print_flows=True, + auxiliary=["ANGLDEGX", "CDIST"], + filename="exchng.gwfgwf", + dev_interfacemodel_on=True, + ) + + # Instantiating MODFLOW 6 solver for flow model + imsgwf = flopy.mf6.ModflowIms( + sim, + print_option="SUMMARY", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="NONE", + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=rclose, + linear_acceleration="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relax, + filename=f"{gwf1.name}.ims", + ) + sim.register_ims_package(imsgwf, [gwf1.name, gwf2.name]) + + # left gwe model + gwe1 = build_gwe_model(idx, sim, "gweleft", side="left") + + # right gwe model + gwe2 = build_gwe_model(idx, sim, "gweright", side="right") + + imsgwe = flopy.mf6.ModflowIms( + sim, + print_option="SUMMARY", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="NONE", + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=rclose, + linear_acceleration="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relax, + filename=f"{gwe1.name}.ims", + ) + sim.register_ims_package(imsgwe, [gwe1.name, gwe2.name]) + + # Instantiating MODFLOW 6 flow-transport exchange mechanism + flopy.mf6.ModflowGwfgwe( + sim, + exgtype="GWF6-GWE6", + exgmnamea="gwfleft", + exgmnameb="gweleft", + pname="GWFGWEl", + filename=f"{gwe1.name}.gwfgwe", + ) + + # Instantiating MODFLOW 6 flow-transport exchange mechanism + flopy.mf6.ModflowGwfgwe( + sim, + exgtype="GWF6-GWE6", + exgmnamea="gwfright", + exgmnameb="gweright", + pname="GWFGWEr", + filename=f"{gwe2.name}.gwfgwe", + ) + + # Create GWE GWE exchange + flopy.mf6.ModflowGwegwe( + sim, + exgtype="GWE6-GWE6", + gwfmodelname1=gwf1.name, + gwfmodelname2=gwf2.name, + adv_scheme=scheme, + nexg=len(exgdata), + exgmnamea=gwe1.name, + exgmnameb=gwe2.name, + exchangedata=exgdata, + auxiliary=["ANGLDEGX", "CDIST"], + filename="exchng.gwegwe", + ) + + return sim, None + + +def check_output(idx, test): + print("evaluating results...") + + # Look for successful termination in simulation lst file + srch_str = "Normal termination of simulation" + fpth = os.path.join(test.workspace, "mfsim.lst") + proceed = False + with open(fpth) as f: + lines = f.readlines() + if any(srch_str in line for line in lines): + proceed = True + + if proceed: + # read transport results from GWE model + gwename1 = "gweleft" + ".ucn" + gwename2 = "gweright" + ".ucn" + + fpth1 = os.path.join(test.workspace, gwename1) + fpth2 = os.path.join(test.workspace, gwename2) + + # load temperatures + tobj1 = flopy.utils.HeadFile( + fpth1, precision="double", text="TEMPERATURE" + ) + tobj2 = flopy.utils.HeadFile( + fpth2, precision="double", text="TEMPERATURE" + ) + temps1 = tobj1.get_alldata() + temps2 = tobj2.get_alldata() + temps_all = np.concatenate((temps1, temps2), axis=3) + np.all(temps_all == strt_temp) + msg0 = "Model finished, but output not uniformly 10.0" + assert np.allclose(temps_all, 10.0), msg0 + else: + msg1 = "Model expected to fail when idx > 0" + print(msg1) + + +# - No need to change any code below +@pytest.mark.parametrize( + "idx, name", + list(enumerate(cases)), +) +def test_mf6model(idx, name, function_tmpdir, targets): + xfail = ["bad" in cases[ct] for ct in np.arange(len(cases))] + test = TestFramework( + name=name, + workspace=function_tmpdir, + targets=targets, + build=lambda t: build_models(idx, t), + check=lambda t: check_output(idx, t), + xfail=xfail[idx], + ) + test.run() From 1e638468a21d307cce821f7c9f3b90170de362f5 Mon Sep 17 00:00:00 2001 From: emorway-usgs Date: Wed, 28 Aug 2024 14:38:47 -0700 Subject: [PATCH 6/7] minor fix --- src/Model/Connection/GweGweConnection.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Model/Connection/GweGweConnection.f90 b/src/Model/Connection/GweGweConnection.f90 index 1edfc2aba92..e4686dda339 100644 --- a/src/Model/Connection/GweGweConnection.f90 +++ b/src/Model/Connection/GweGweConnection.f90 @@ -492,7 +492,7 @@ end subroutine gwegwecon_ot !< subroutine validateGweExchange(this) ! modules - use SimVariablesModule, only: errmsg, simulation_mode + use SimVariablesModule, only: errmsg use SimModule, only: store_error use GweEstModule, only: GweEstType From b0b335b9615cb7b962e5d4e3c330607348cb0285 Mon Sep 17 00:00:00 2001 From: emorway-usgs Date: Wed, 28 Aug 2024 15:03:12 -0700 Subject: [PATCH 7/7] always something --- src/Model/Connection/GweGweConnection.f90 | 112 +++++++++++----------- 1 file changed, 56 insertions(+), 56 deletions(-) diff --git a/src/Model/Connection/GweGweConnection.f90 b/src/Model/Connection/GweGweConnection.f90 index e4686dda339..b7cd571e282 100644 --- a/src/Model/Connection/GweGweConnection.f90 +++ b/src/Model/Connection/GweGweConnection.f90 @@ -93,24 +93,24 @@ subroutine gweGweConnection_ctor(this, model, gweEx) character(len=LENCOMPONENTNAME) :: name class(*), pointer :: objPtr logical(LGP) :: write_ifmodel_listfile = .false. - + objPtr => model this%gweModel => CastAsGweModel(objPtr) objPtr => gweEx this%gweExchange => CastAsGweExchange(objPtr) - + if (gweEx%v_model1%is_local .and. gweEx%v_model2%is_local) then this%owns_exchange = associated(model, gweEx%model1) else this%owns_exchange = .true. end if - + if (gweEx%v_model1 == model) then write (name, '(a,i0)') 'GWECON1_', gweEx%id else write (name, '(a,i0)') 'GWECON2_', gweEx%id end if - + ! .lst file for interface model if (write_ifmodel_listfile) then fname = trim(name)//'.im.lst' @@ -119,21 +119,21 @@ subroutine gweGweConnection_ctor(this, model, gweEx) trim(this%gweModel%name), 'from exchange ', & trim(gweEx%name) end if - + ! first call base constructor call this%SpatialModelConnectionType%spatialConnection_ctor(model, & gweEx, & name) - + call this%allocate_scalars() this%typename = 'GWE-GWE' this%iIfaceAdvScheme = 0 this%iIfaceXt3d = 0 this%exgflowSign = 1 - + allocate (this%gweInterfaceModel) this%interface_model => this%gweInterfaceModel - + end subroutine gweGweConnection_ctor !> @brief Allocate scalar variables for this connection @@ -141,11 +141,11 @@ end subroutine gweGweConnection_ctor subroutine allocate_scalars(this) ! dummy class(GweGweConnectionType) :: this !< the connection - + call mem_allocate(this%iIfaceAdvScheme, 'IADVSCHEME', this%memoryPath) call mem_allocate(this%iIfaceXt3d, 'IXT3D', this%memoryPath) call mem_allocate(this%exgflowSign, 'EXGFLOWSIGN', this%memoryPath) - + end subroutine allocate_scalars !> @brief define the GWE-GWE connection @@ -159,7 +159,7 @@ subroutine gwegwecon_df(this) ! determine advection scheme (the GWE-GWE exchange ! has been read at this point) this%iIfaceAdvScheme = this%gweExchange%iAdvScheme - + ! determine xt3d setting on interface this%iIfaceXt3d = this%gweExchange%ixt3d @@ -211,14 +211,14 @@ end subroutine gwegwecon_df subroutine cfg_dist_vars(this) ! dummy class(GweGweConnectionType) :: this !< the connection - + call this%cfg_dv('X', '', SYNC_NDS, & (/STG_BFR_CON_AR, STG_BFR_EXG_AD, STG_BFR_EXG_CF/)) call this%cfg_dv('IBOUND', '', SYNC_NDS, (/STG_BFR_CON_AR/)) call this%cfg_dv('TOP', 'DIS', SYNC_NDS, (/STG_BFR_CON_AR/)) call this%cfg_dv('BOT', 'DIS', SYNC_NDS, (/STG_BFR_CON_AR/)) call this%cfg_dv('AREA', 'DIS', SYNC_NDS, (/STG_BFR_CON_AR/)) - + if (this%gweInterfaceModel%cnd%idisp > 0) then call this%cfg_dv('ALH', 'CND', SYNC_NDS, (/STG_BFR_CON_AR/)) call this%cfg_dv('ALV', 'CND', SYNC_NDS, (/STG_BFR_CON_AR/)) @@ -238,7 +238,7 @@ subroutine cfg_dist_vars(this) if (this%gweModel%incnd > 0 .and. this%gweModel%inest > 0) then call this%cfg_dv('POROSITY', 'EST', SYNC_NDS, (/STG_AFT_CON_AR/)) end if - + end subroutine cfg_dist_vars !> @brief Allocate array variables for this connection @@ -258,10 +258,10 @@ subroutine setGridExtent(this) class(GweGweConnectionType) :: this !< the connection ! local logical(LGP) :: hasAdv, hasCnd - + hasAdv = this%gweModel%inadv > 0 hasCnd = this%gweModel%incnd > 0 - + if (hasAdv) then if (this%iIfaceAdvScheme == 2) then this%exg_stencil_depth = 2 @@ -270,7 +270,7 @@ subroutine setGridExtent(this) end if end if end if - + if (hasCnd) then if (this%iIfaceXt3d > 0) then this%exg_stencil_depth = 2 @@ -279,7 +279,7 @@ subroutine setGridExtent(this) end if end if end if - + end subroutine setGridExtent !> @brief allocate and read/set the connection's data structures @@ -336,7 +336,7 @@ subroutine validateConnection(this) ! we cannot validate this (yet) in parallel mode if (.not. this%gweExchange%v_model1%is_local) return if (.not. this%gweExchange%v_model2%is_local) return - + ! check specific cross-interface options/values that should be the same call this%validateGweExchange() @@ -350,7 +350,7 @@ subroutine validateConnection(this) &and the other one is not' call store_error(errmsg) end if - + if ((this%gweExchange%gwemodel1%incnd > 0 .and. & this%gweExchange%gwemodel2%incnd == 0) .or. & (this%gweExchange%gwemodel2%incnd > 0 .and. & @@ -360,39 +360,39 @@ subroutine validateConnection(this) &and the other one is not' call store_error(errmsg) end if - + ! abort on errors if (count_errors() > 0) then write (errmsg, '(a)') 'Errors occurred while processing exchange(s)' call ustop() end if - + end subroutine validateConnection subroutine gwegwecon_rp(this) ! dummy class(GweGweConnectionType) :: this !< the connection - + ! call exchange rp routines if (this%owns_exchange) then call this%gweExchange%exg_rp() end if - + end subroutine gwegwecon_rp !> @brief Advance this connection !< subroutine gwegwecon_ad(this) - + class(GweGweConnectionType) :: this !< this connection - + ! recalculate conduction ellipse if (this%gweInterfaceModel%incnd > 0) call this%gweInterfaceModel%cnd%cnd_ad() - + if (this%owns_exchange) then call this%gweExchange%exg_ad() end if - + end subroutine gwegwecon_ad subroutine gwegwecon_fc(this, kiter, matrix_sln, rhs_sln, inwtflag) @@ -402,10 +402,10 @@ subroutine gwegwecon_fc(this, kiter, matrix_sln, rhs_sln, inwtflag) class(MatrixBaseType), pointer :: matrix_sln !< the system matrix real(DP), dimension(:), intent(inout) :: rhs_sln !< global right-hand-side integer(I4B), optional, intent(in) :: inwtflag !< newton-raphson flag - + call this%SpatialModelConnectionType%spatialcon_fc( & kiter, matrix_sln, rhs_sln, inwtflag) - + ! _fc the movers through the exchange if (this%owns_exchange) then if (this%gweExchange%inmvt > 0) then @@ -413,7 +413,7 @@ subroutine gwegwecon_fc(this, kiter, matrix_sln, rhs_sln, inwtflag) this%gweExchange%gwemodel2%x) end if end if - + end subroutine gwegwecon_fc subroutine gwegwecon_cq(this, icnvg, isuppress_output, isolnid) @@ -422,10 +422,10 @@ subroutine gwegwecon_cq(this, icnvg, isuppress_output, isolnid) integer(I4B), intent(inout) :: icnvg !< convergence flag integer(I4B), intent(in) :: isuppress_output !< suppress output when =1 integer(I4B), intent(in) :: isolnid !< solution id - + call this%gweInterfaceModel%model_cq(icnvg, isuppress_output) call this%setFlowToExchange() - + end subroutine gwegwecon_cq !> @brief Set the flows (flowja from interface model) to the @@ -439,11 +439,11 @@ subroutine setFlowToExchange(this) integer(I4B) :: i class(GweExchangeType), pointer :: gweEx type(IndexMapSgnType), pointer :: map - + if (this%owns_exchange) then gweEx => this%gweExchange map => this%interface_map%exchange_maps(this%interface_map%prim_exg_idx) - + ! use (half of) the exchange map in reverse: do i = 1, size(map%src_idx) if (map%sign(i) < 0) cycle ! simvals is defined from exg%m1 => exg%m2 @@ -451,7 +451,7 @@ subroutine setFlowToExchange(this) this%gweInterfaceModel%flowja(map%tgt_idx(i)) end do end if - + end subroutine setFlowToExchange subroutine gwegwecon_bd(this, icnvg, isuppress_output, isolnid) @@ -462,26 +462,26 @@ subroutine gwegwecon_bd(this, icnvg, isuppress_output, isolnid) integer(I4B), intent(inout) :: icnvg !< convergence flag integer(I4B), intent(in) :: isuppress_output !< suppress output when =1 integer(I4B), intent(in) :: isolnid !< solution id - + ! call exchange budget routine, also calls _bd ! for movers. if (this%owns_exchange) then call this%gweExchange%exg_bd(icnvg, isuppress_output, isolnid) end if - + end subroutine gwegwecon_bd subroutine gwegwecon_ot(this) ! dummy class(GweGweConnectionType) :: this !< the connection - + ! call exg_ot() here as it handles all output processing ! based on gweExchange%simvals(:), which was correctly ! filled from gwegwecon if (this%owns_exchange) then call this%gweExchange%exg_ot() end if - + end subroutine gwegwecon_ot !> @brief Validate the exchange, intercepting those @@ -495,10 +495,10 @@ subroutine validateGweExchange(this) use SimVariablesModule, only: errmsg use SimModule, only: store_error use GweEstModule, only: GweEstType - + ! dummy class(GweGweConnectionType) :: this !< this connection - + ! local class(GweExchangeType), pointer :: gweEx class(*), pointer :: modelPtr @@ -517,7 +517,7 @@ subroutine validateGweExchange(this) gweModel1 => CastAsGweModel(modelPtr) modelPtr => this%gweExchange%model2 gweModel2 => CastAsGweModel(modelPtr) - + ! check that EST package usage is the same on both side of the interface if ((gweModel1%inest > 0 .and. gweModel2%inest == 0) .or. & (gweModel1%inest == 0 .and. gweModel2%inest > 0)) then @@ -526,14 +526,14 @@ subroutine validateGweExchange(this) 'interface model for exchange ', & trim(gweEx%name) call store_error(errmsg) - end if + end if ! conduction options need to be the same in both model if ((gweModel1%cnd%ixt3d > 0 .and. gweModel2%cnd%ixt3d == 0) .or. & (gweModel1%cnd%ixt3d == 0 .and. gweModel2%cnd%ixt3d > 0)) then write (errmsg, '(2a)') 'Use of XT3D to calculate conduction should '// & 'be the same in both models, either both use XT3D or neither for '// & - ' exchange '// trim(gweEx%name) + ' exchange '//trim(gweEx%name) call store_error(errmsg) end if @@ -551,7 +551,7 @@ subroutine validateGweExchange(this) !end if if (.not. compatible) then write (errmsg, '(6a)') 'Energy storage and transfer (EST) packages ', & - 'in model '// trim(gweEx%model1%name), ' and ', & + 'in model '//trim(gweEx%model1%name), ' and ', & trim(gweEx%model2%name), & ' should be equivalent to construct an '// & ' interface model for exchange ', & @@ -568,32 +568,32 @@ subroutine gwegwecon_da(this) class(GweGweConnectionType) :: this !< the connection ! local logical(LGP) :: isOpen - + ! scalars call mem_deallocate(this%iIfaceAdvScheme) call mem_deallocate(this%iIfaceXt3d) call mem_deallocate(this%exgflowSign) - + ! arrays call mem_deallocate(this%exgflowjaGwe) - + ! interface model call this%gweInterfaceModel%model_da() deallocate (this%gweInterfaceModel) - + ! dealloc base call this%spatialcon_da() - + inquire (this%iout, opened=isOpen) if (isOpen) then close (this%iout) end if - + ! we need to deallocate the exchange we own: if (this%owns_exchange) then call this%gweExchange%exg_da() end if - + end subroutine gwegwecon_da !> @brief Cast to GweGweConnectionType @@ -604,15 +604,15 @@ function CastAsGweGweConnection(obj) result(res) class(*), pointer, intent(inout) :: obj !< object to be cast ! return class(GweGweConnectionType), pointer :: res !< the GweGweConnection - + res => null() if (.not. associated(obj)) return - + select type (obj) class is (GweGweConnectionType) res => obj end select - + end function CastAsGweGweConnection end module