From d41878c492f0d1c21fa4d610eccd2f04a2dcd9aa Mon Sep 17 00:00:00 2001 From: Eric Morway Date: Thu, 12 Oct 2023 07:50:33 -0700 Subject: [PATCH] refactor(tsp): Elevate IC to generalized transport class (#1392) --- make/makefile | 4 +- msvs/mf6core.vfproj | 4 +- src/Model/GroundWaterTransport/gwt1.f90 | 85 +++++++++++-------- src/Model/GroundWaterTransport/gwt1obs1.f90 | 8 +- src/Model/TransportModel/tsp1.f90 | 16 +++- .../tsp1ic1.f90} | 50 +++++------ src/meson.build | 2 +- 7 files changed, 99 insertions(+), 70 deletions(-) rename src/Model/{GroundWaterTransport/gwt1ic1.f90 => TransportModel/tsp1ic1.f90} (68%) diff --git a/make/makefile b/make/makefile index 0afe8256a21..6867c3485ae 100644 --- a/make/makefile +++ b/make/makefile @@ -206,11 +206,11 @@ $(OBJDIR)/gwf3drn8.o \ $(OBJDIR)/IndexMap.o \ $(OBJDIR)/VirtualModel.o \ $(OBJDIR)/BaseExchange.o \ +$(OBJDIR)/gwf3ic8.o \ $(OBJDIR)/tsp1fmi1.o \ $(OBJDIR)/TspAdvOptions.o \ $(OBJDIR)/UzfCellGroup.o \ $(OBJDIR)/OutputControlData.o \ -$(OBJDIR)/gwf3ic8.o \ $(OBJDIR)/Xt3dInterface.o \ $(OBJDIR)/gwf3tvk8.o \ $(OBJDIR)/gwf3vsc8.o \ @@ -221,6 +221,7 @@ $(OBJDIR)/ImsLinearSettings.o \ $(OBJDIR)/ConvergenceSummary.o \ $(OBJDIR)/CellWithNbrs.o \ $(OBJDIR)/NumericalExchange.o \ +$(OBJDIR)/tsp1ic1.o \ $(OBJDIR)/tsp1adv1.o \ $(OBJDIR)/gwf3disv8.o \ $(OBJDIR)/gwf3disu8.o \ @@ -229,7 +230,6 @@ $(OBJDIR)/gwf3uzf8.o \ $(OBJDIR)/gwt1apt1.o \ $(OBJDIR)/GwtSpc.o \ $(OBJDIR)/OutputControl.o \ -$(OBJDIR)/gwt1ic1.o \ $(OBJDIR)/gwt1mst1.o \ $(OBJDIR)/GwtDspOptions.o \ $(OBJDIR)/gwf3npf8.o \ diff --git a/msvs/mf6core.vfproj b/msvs/mf6core.vfproj index 55ea9d954e8..e9cd6a8356d 100644 --- a/msvs/mf6core.vfproj +++ b/msvs/mf6core.vfproj @@ -170,7 +170,6 @@ - @@ -208,7 +207,8 @@ - + + diff --git a/src/Model/GroundWaterTransport/gwt1.f90 b/src/Model/GroundWaterTransport/gwt1.f90 index 901ab373f1d..50c5e8c5941 100644 --- a/src/Model/GroundWaterTransport/gwt1.f90 +++ b/src/Model/GroundWaterTransport/gwt1.f90 @@ -15,7 +15,6 @@ module GwtModule use TransportModelModule, only: TransportModelType use BaseModelModule, only: BaseModelType use BndModule, only: BndType, AddBndToList, GetBndFromList - use GwtIcModule, only: GwtIcType use TspFmiModule, only: TspFmiType use GwtDspModule, only: GwtDspType use GwtSsmModule, only: GwtSsmType @@ -38,14 +37,12 @@ module GwtModule type, extends(TransportModelType) :: GwtModelType - type(GwtIcType), pointer :: ic => null() ! initial conditions package type(GwtMstType), pointer :: mst => null() ! mass storage and transfer package type(GwtDspType), pointer :: dsp => null() ! dispersion package type(GwtSsmType), pointer :: ssm => null() ! source sink mixing package type(GwtMvtType), pointer :: mvt => null() ! mover transport package type(GwtOcType), pointer :: oc => null() ! output control package type(GwtObsType), pointer :: obs => null() ! observation package - integer(I4B), pointer :: inic => null() ! unit number IC integer(I4B), pointer :: inmvt => null() ! unit number MVT integer(I4B), pointer :: inmst => null() ! unit number MST integer(I4B), pointer :: indsp => null() ! DSP enabled flag @@ -158,10 +155,11 @@ subroutine gwt_cr(filename, id, modelname) return end subroutine gwt_cr - !> @brief Define packages of the model - ! - ! (1) call df routines for each package - ! (2) set variables and pointers + !> @brief Define packages of the GWT model + !! + !! This subroutine defines a gwt model type. Steps include: + !! (1) call df routines for each package + !! (2) set variables and pointers !< subroutine gwt_df(this) ! -- modules @@ -243,7 +241,8 @@ subroutine gwt_ac(this, sparse) return end subroutine gwt_ac - !> @brief Map connection positions in numerical solution coefficient matrix. + !> @brief Map the positions of the GWT model connections in the numerical + !! solution coefficient matrix. !< subroutine gwt_mc(this, matrix_sln) ! -- dummy @@ -256,6 +255,7 @@ subroutine gwt_mc(this, matrix_sln) ! -- Find the position of each connection in the global ia, ja structure ! and store them in idxglo. call this%dis%dis_mc(this%moffset, this%idxglo, matrix_sln) + ! if (this%indsp > 0) call this%dsp%dsp_mc(this%moffset, matrix_sln) ! ! -- Map any package connections @@ -268,10 +268,11 @@ subroutine gwt_mc(this, matrix_sln) return end subroutine gwt_mc - !> @brief Allocate and Read - ! - ! (1) allocates and reads packages part of this model, - ! (2) allocates memory for arrays part of this model object + !> @brief GWT Model Allocate and Read + !! + !! This subroutine: + !! - allocates and reads packages that are part of this model, + !! - allocates memory for arrays used by this model object !< subroutine gwt_ar(this) ! -- modules @@ -321,7 +322,9 @@ subroutine gwt_ar(this) return end subroutine gwt_ar - !> @brief Read and prepare (calls package read and prepare routines) + !> @brief GWT Model Read and Prepare + !! + !! Call the read and prepare routines of the attached packages !< subroutine gwt_rp(this) ! -- modules @@ -352,7 +355,9 @@ subroutine gwt_rp(this) return end subroutine gwt_rp - !> @brief Time step advance (calls package advance subroutines) + !> @brief GWT Model Time Step Advance + !! + !! Call the advance subroutines of the attached packages !< subroutine gwt_ad(this) ! -- modules @@ -407,9 +412,12 @@ subroutine gwt_ad(this) return end subroutine gwt_ad - !> @brief Calculate coefficients + !> @brief GWT Model calculate coefficients + !! + !! Call the calculate coefficients subroutines of the attached packages !< subroutine gwt_cf(this, kiter) + ! -- modules ! -- dummy class(GwtModelType) :: this integer(I4B), intent(in) :: kiter @@ -427,9 +435,12 @@ subroutine gwt_cf(this, kiter) return end subroutine gwt_cf - !> @brief Fill coefficients + !> @brief GWT Model fill coefficients + !! + !! Call the fill coefficients subroutines attached packages !< subroutine gwt_fc(this, kiter, matrix_sln, inwtflag) + ! -- modules ! -- dummy class(GwtModelType) :: this integer(I4B), intent(in) :: kiter @@ -471,7 +482,9 @@ subroutine gwt_fc(this, kiter, matrix_sln, inwtflag) return end subroutine gwt_fc - !> @brief Final convergence check (calls package cc routines) + !> @brief GWT Model Final Convergence Check + !! + !! If MVR/MVT is active, call the MVR convergence check subroutines !< subroutine gwt_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak) ! -- dummy @@ -501,7 +514,9 @@ subroutine gwt_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak) return end subroutine gwt_cc - !> @brief Calculate intercell flows (flowja) + !> @brief GWT Model calculate flow + !! + !! Call the intercell flows (flow ja) subroutine !< subroutine gwt_cq(this, icnvg, isuppress_output) ! -- modules @@ -548,10 +563,11 @@ subroutine gwt_cq(this, icnvg, isuppress_output) return end subroutine gwt_cq - !> @brief Model budget - ! - ! (1) Calculate intercell flows (flowja) - ! (2) Calculate package contributions to model budget + !> @brief GWT Model Budget + !! + !! This subroutine: + !! (1) calculates intercell flows (flowja) + !! (2) calculates package contributions to the model budget !< subroutine gwt_bd(this, icnvg, isuppress_output) use ConstantsModule, only: DZERO @@ -579,13 +595,14 @@ subroutine gwt_bd(this, icnvg, isuppress_output) packobj => GetBndFromList(this%bndlist, ip) call packobj%bnd_bd(this%budget) end do - ! ! -- Return return end subroutine gwt_bd !> @brief Print and/or save model output + !! + !! Call the parent class output routine !< subroutine gwt_ot(this) ! -- modules @@ -813,6 +830,8 @@ subroutine gwt_ot_bdsummary(this, ibudfl, ipflag) end subroutine gwt_ot_bdsummary !> @brief Deallocate + !! + !! Deallocate memmory at conclusion of model run !< subroutine gwt_da(this) ! -- modules @@ -862,7 +881,6 @@ subroutine gwt_da(this) end do ! ! -- Scalars - call mem_deallocate(this%inic) call mem_deallocate(this%indsp) call mem_deallocate(this%inssm) call mem_deallocate(this%inmst) @@ -885,8 +903,6 @@ end subroutine gwt_da !! This subroutine adds a budget entry to the flow budget. It was added as !! a method for the gwt model object so that the exchange object could add its !! contributions. - !! - !! (1) adds the entry to the budget object !< subroutine gwt_bdentry(this, budterm, budtxt, rowlabel) ! -- modules @@ -937,7 +953,11 @@ function gwt_get_iasym(this) result(iasym) return end function gwt_get_iasym - !> @brief Allocate memory for non-allocatable members + !> Allocate memory for non-allocatable members + !! + !! A subroutine for allocating the scalars specific to the GWT model type. + !! Additional scalars used by the parent class are allocated by the parent + !! class. !< subroutine allocate_scalars(this, modelname) ! -- modules @@ -950,7 +970,6 @@ subroutine allocate_scalars(this, modelname) call this%allocate_tsp_scalars(modelname) ! ! -- allocate members that are part of model class - call mem_allocate(this%inic, 'INIC', this%memoryPath) call mem_allocate(this%inmvt, 'INMVT', this%memoryPath) call mem_allocate(this%inmst, 'INMST', this%memoryPath) call mem_allocate(this%indsp, 'INDSP', this%memoryPath) @@ -958,7 +977,6 @@ subroutine allocate_scalars(this, modelname) call mem_allocate(this%inoc, 'INOC ', this%memoryPath) call mem_allocate(this%inobs, 'INOBS', this%memoryPath) ! - this%inic = 0 this%inmvt = 0 this%inmst = 0 this%indsp = 0 @@ -971,6 +989,8 @@ subroutine allocate_scalars(this, modelname) end subroutine allocate_scalars !> @brief Create boundary condition packages for this model + !! + !! Call the package create routines for packages activated by the user. !< subroutine package_create(this, filtyp, ipakid, ipaknum, pakname, inunit, & iout) @@ -1056,6 +1076,7 @@ function CastAsGwtModel(model) result(gwtmodel) type is (GwtModelType) gwtmodel => model end select + ! ! -- Return return end function CastAsGwtModel @@ -1127,7 +1148,6 @@ subroutine create_gwt_packages(this, indis) use MemoryManagerModule, only: mem_setptr use MemoryHelperModule, only: create_mem_path use SimVariablesModule, only: idm_context - use GwtIcModule, only: ic_cr use GwtMstModule, only: mst_cr use GwtDspModule, only: dsp_cr use GwtSsmModule, only: ssm_cr @@ -1174,8 +1194,6 @@ subroutine create_gwt_packages(this, indis) ! ! -- create dis package first as it is a prerequisite for other packages select case (pkgtype) - case ('IC6') - this%inic = inunit case ('MVT6') this%inmvt = inunit case ('MST6') @@ -1199,7 +1217,6 @@ subroutine create_gwt_packages(this, indis) end do ! ! -- Create packages that are tied directly to model - call ic_cr(this%ic, this%name, this%inic, this%iout, this%dis) call mst_cr(this%mst, this%name, this%inmst, this%iout, this%fmi) call dsp_cr(this%dsp, this%name, mempathdsp, this%indsp, this%iout, & this%fmi) @@ -1209,7 +1226,7 @@ subroutine create_gwt_packages(this, indis) call gwt_obs_cr(this%obs, this%inobs) ! ! -- Check to make sure that required ftype's have been specified - call this%ftype_check(indis, this%inmst, this%inic) + call this%ftype_check(indis, this%inmst) ! call this%create_bndpkgs(bndpkgs, pkgtypes, pkgnames, mempaths, inunits) ! diff --git a/src/Model/GroundWaterTransport/gwt1obs1.f90 b/src/Model/GroundWaterTransport/gwt1obs1.f90 index 48dd58f0e7c..fc81782f491 100644 --- a/src/Model/GroundWaterTransport/gwt1obs1.f90 +++ b/src/Model/GroundWaterTransport/gwt1obs1.f90 @@ -3,7 +3,7 @@ module GwtObsModule use KindModule, only: DP, I4B use ConstantsModule, only: LINELENGTH, MAXOBSTYPES use BaseDisModule, only: DisBaseType - use GwtIcModule, only: GwtIcType + use TspIcModule, only: TspIcType use ObserveModule, only: ObserveType use ObsModule, only: ObsType use SimModule, only: count_errors, store_error, & @@ -15,7 +15,7 @@ module GwtObsModule type, extends(ObsType) :: GwtObsType ! -- Private members - type(GwtIcType), pointer, private :: ic => null() ! initial conditions + type(TspIcType), pointer, private :: ic => null() ! initial conditions real(DP), dimension(:), pointer, contiguous, private :: x => null() ! concentration real(DP), dimension(:), pointer, contiguous, private :: flowja => null() ! intercell flows contains @@ -64,7 +64,7 @@ subroutine gwt_obs_ar(this, ic, x, flowja) ! ------------------------------------------------------------------------------ ! -- dummy class(GwtObsType), intent(inout) :: this - type(GwtIcType), pointer, intent(in) :: ic + type(TspIcType), pointer, intent(in) :: ic real(DP), dimension(:), pointer, contiguous, intent(in) :: x real(DP), dimension(:), pointer, contiguous, intent(in) :: flowja ! ------------------------------------------------------------------------------ @@ -193,7 +193,7 @@ subroutine set_pointers(this, ic, x, flowja) ! ------------------------------------------------------------------------------ ! -- dummy class(GwtObsType), intent(inout) :: this - type(GwtIcType), pointer, intent(in) :: ic + type(TspIcType), pointer, intent(in) :: ic real(DP), dimension(:), pointer, contiguous, intent(in) :: x real(DP), dimension(:), pointer, contiguous, intent(in) :: flowja ! ------------------------------------------------------------------------------ diff --git a/src/Model/TransportModel/tsp1.f90 b/src/Model/TransportModel/tsp1.f90 index 124036429c8..239f3a2a5d1 100644 --- a/src/Model/TransportModel/tsp1.f90 +++ b/src/Model/TransportModel/tsp1.f90 @@ -11,6 +11,7 @@ module TransportModelModule LENMEMPATH, LENVARNAME use SimVariablesModule, only: errmsg use NumericalModelModule, only: NumericalModelType + use TspIcModule, only: TspIcType use TspFmiModule, only: TspFmiType use TspAdvModule, only: TspAdvType use BudgetModule, only: BudgetType @@ -26,9 +27,11 @@ module TransportModelModule type(TspFmiType), pointer :: fmi => null() ! flow model interface type(TspAdvType), pointer :: adv => null() !< advection package + type(TspIcType), pointer :: ic => null() !< initial conditions package type(BudgetType), pointer :: budget => null() !< budget object integer(I4B), pointer :: infmi => null() ! unit number FMI integer(I4B), pointer :: inadv => null() !< unit number ADV + integer(I4B), pointer :: inic => null() !< unit number IC real(DP), pointer :: eqnsclfac => null() !< constant factor by which all terms in the model's governing equation are scaled (divided) for formulation and solution ! Labels that will be defined character(len=LENVARNAME) :: tsptype = '' !< "solute" or "heat" @@ -279,10 +282,12 @@ subroutine allocate_tsp_scalars(this, modelname) call this%NumericalModelType%allocate_scalars(modelname) ! ! -- allocate members that are part of model class + call mem_allocate(this%inic, 'INIC', this%memoryPath) call mem_allocate(this%infmi, 'INFMI', this%memoryPath) call mem_allocate(this%inadv, 'INADV', this%memoryPath) call mem_allocate(this%eqnsclfac, 'EQNSCLFAC', this%memoryPath) ! + this%inic = 0 this%infmi = 0 this%inadv = 0 this%eqnsclfac = DZERO @@ -332,6 +337,7 @@ subroutine tsp_da(this) ! -- local ! ! -- Scalars + call mem_deallocate(this%inic) call mem_deallocate(this%infmi) call mem_deallocate(this%inadv) call mem_deallocate(this%eqnsclfac) @@ -344,7 +350,7 @@ end subroutine tsp_da !! !! Check to make sure required input files have been specified !< - subroutine ftype_check(this, indis, inmst, inic) + subroutine ftype_check(this, indis, inmst) ! -- modules use ConstantsModule, only: LINELENGTH use SimModule, only: store_error, count_errors, store_error_filename @@ -352,13 +358,12 @@ subroutine ftype_check(this, indis, inmst, inic) class(TransportModelType) :: this integer(I4B), intent(in) :: indis integer(I4B), intent(in) :: inmst - integer(I4B), intent(in) :: inic ! -- local character(len=LINELENGTH) :: errmsg ! ------------------------------------------------------------------------------ ! ! -- Check for IC6, DIS(u), and MST. Stop if not present. - if (inic == 0) then + if (this%inic == 0) then write (errmsg, '(a)') & 'Initial conditions (IC6) package not specified.' call store_error(errmsg) @@ -487,6 +492,7 @@ subroutine create_tsp_packages(this, indis) use GwfDisModule, only: dis_cr use GwfDisvModule, only: disv_cr use GwfDisuModule, only: disu_cr + use TspIcModule, only: ic_cr use TspFmiModule, only: fmi_cr use TspAdvModule, only: adv_cr ! -- dummy @@ -539,6 +545,8 @@ subroutine create_tsp_packages(this, indis) case ('DISU6') indis = 1 call disu_cr(this%dis, this%name, mempath, indis, this%iout) + case ('IC6') + this%inic = inunit case ('FMI6') this%infmi = inunit case ('ADV6') @@ -547,6 +555,8 @@ subroutine create_tsp_packages(this, indis) end do ! ! -- Create packages that are tied directly to model + call ic_cr(this%ic, this%name, this%inic, this%iout, this%dis, & + this%depvartype) call fmi_cr(this%fmi, this%name, this%infmi, this%iout, this%eqnsclfac, & this%depvartype) call adv_cr(this%adv, this%name, this%inadv, this%iout, this%fmi, & diff --git a/src/Model/GroundWaterTransport/gwt1ic1.f90 b/src/Model/TransportModel/tsp1ic1.f90 similarity index 68% rename from src/Model/GroundWaterTransport/gwt1ic1.f90 rename to src/Model/TransportModel/tsp1ic1.f90 index e9d872a7137..cb060cfb658 100644 --- a/src/Model/GroundWaterTransport/gwt1ic1.f90 +++ b/src/Model/TransportModel/tsp1ic1.f90 @@ -1,37 +1,39 @@ -module GwtIcModule +module TspIcModule use KindModule, only: DP, I4B + use ConstantsModule, only: LENVARNAME use GwfIcModule, only: GwfIcType use BlockParserModule, only: BlockParserType use BaseDisModule, only: DisBaseType implicit none private - public :: GwtIcType + public :: TspIcType public :: ic_cr - ! -- Most of the GwtIcType functionality comes from GwfIcType - type, extends(GwfIcType) :: GwtIcType + ! -- Most of the TspIcType functionality comes from GwfIcType + type, extends(GwfIcType) :: TspIcType + ! -- strings + character(len=LENVARNAME) :: depvartype = '' + contains + procedure :: read_data - end type GwtIcType + + end type TspIcType contains - subroutine ic_cr(ic, name_model, inunit, iout, dis) -! ****************************************************************************** -! ic_cr -- Create a new initial conditions object -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ + !> @brief Create a new initial conditions object + !< + subroutine ic_cr(ic, name_model, inunit, iout, dis, depvartype) ! -- dummy - type(GwtIcType), pointer :: ic + type(TspIcType), pointer :: ic character(len=*), intent(in) :: name_model integer(I4B), intent(in) :: inunit integer(I4B), intent(in) :: iout class(DisBaseType), pointer, intent(in) :: dis -! ------------------------------------------------------------------------------ + character(len=LENVARNAME), intent(in) :: depvartype ! ! -- Create the object allocate (ic) @@ -48,6 +50,9 @@ subroutine ic_cr(ic, name_model, inunit, iout, dis) ! -- set pointers ic%dis => dis ! + ! -- Give package access to the assigned labelsd based on dependent variable + ic%depvartype = depvartype + ! ! -- Initialize block parser call ic%parser%Initialize(ic%inunit, ic%iout) ! @@ -55,18 +60,16 @@ subroutine ic_cr(ic, name_model, inunit, iout, dis) return end subroutine ic_cr + !> @brief Read initial conditions + !! + !! Read initial concentrations or temperatures depending on model type + !< subroutine read_data(this) -! ****************************************************************************** -! read_data -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use ConstantsModule, only: LINELENGTH use SimModule, only: store_error ! -- dummy - class(GwtIcType) :: this + class(TspIcType) :: this ! -- local character(len=LINELENGTH) :: errmsg, keyword character(len=:), allocatable :: line @@ -74,10 +77,9 @@ subroutine read_data(this) logical :: isfound, endOfBlock character(len=24) :: aname(1) ! -- formats -! ------------------------------------------------------------------------------ ! ! -- Setup the label - aname(1) = 'INITIAL CONCENTRATION' + write (aname(1), '(a,1x,a)') 'INITIAL', trim(adjustl(this%depvartype)) ! ! -- get griddata block call this%parser%GetBlock('GRIDDATA', isfound, ierr) @@ -111,4 +113,4 @@ subroutine read_data(this) return end subroutine read_data -end module GwtIcModule +end module TspIcModule diff --git a/src/meson.build b/src/meson.build index 6b84eab5c8f..d72ab2f9b7f 100644 --- a/src/meson.build +++ b/src/meson.build @@ -97,7 +97,6 @@ modflow_sources = files( 'Model' / 'GroundWaterTransport' / 'gwt1disv1idm.f90', 'Model' / 'GroundWaterTransport' / 'gwt1dsp1.f90', 'Model' / 'GroundWaterTransport' / 'gwt1dsp1idm.f90', - 'Model' / 'GroundWaterTransport' / 'gwt1ic1.f90', 'Model' / 'GroundWaterTransport' / 'gwt1idm.f90', 'Model' / 'GroundWaterTransport' / 'gwt1ist1.f90', 'Model' / 'GroundWaterTransport' / 'gwt1lkt1.f90', @@ -134,6 +133,7 @@ modflow_sources = files( 'Model' / 'TransportModel' / 'tsp1.f90', 'Model' / 'TransportModel' / 'tsp1adv1.f90', 'Model' / 'TransportModel' / 'tsp1fmi1.f90', + 'Model' / 'TransportModel' / 'tsp1ic1.f90', 'Model' / 'BaseModel.f90', 'Model' / 'ExplicitModel.f90', 'Model' / 'NumericalModel.f90',