diff --git a/autotest/simulation.py b/autotest/simulation.py index 9ee16586398..81e97eb460c 100644 --- a/autotest/simulation.py +++ b/autotest/simulation.py @@ -17,7 +17,7 @@ from flopy.utils.compare import compare_heads from modflow_devtools.misc import is_in_ci -DNODATA = 3.e+30 +DNODATA = 3.0e30 sfmt = "{:25s} - {}" extdict = { "hds": "head", @@ -164,7 +164,6 @@ def setup(self, src, dst): return def setup_comparison(self, src, dst, testModel=True): - # evaluate if comparison should be made if not self.make_comparison: return @@ -221,7 +220,7 @@ def run(self): success, buff = self.run_parallel( exe, ) - except Exception as exc: + except Exception as exc: msg = sfmt.format("MODFLOW 6 run", self.name) print(msg) print(exc) @@ -332,12 +331,37 @@ def run(self): return def run_parallel(self, exe): - normal_msg="normal termination" + physical_cpus = os.cpu_count() + print(f"CPUs: {physical_cpus}") + if self.ncpus > physical_cpus: + print( + f"simulation is oversubscribed to {self.ncpus} CPUs " + + f"but there are only {physical_cpus} CPUs. " + + "Expect degraded performance." + ) + is_oversubscribed = True + with open(f"{self.simpath}/localhost", "w") as f: + f.write(f"localhost slots={self.ncpus}\n") + else: + is_oversubscribed = False + + normal_msg = "normal termination" success = False nr_success = 0 buff = [] - mpiexec_cmd = ["mpiexec", "--oversubscribe", "-np", str(self.ncpus), exe, "-p"] + # add initial parallel commands + mpiexec_cmd = ["mpiexec", "-np", str(self.ncpus)] + + # add oversubscribed commands + if is_oversubscribed: + mpiexec_cmd.append("--hostfile") + mpiexec_cmd.append("localhost") + + # add remainder of parallel commands + mpiexec_cmd.append(exe) + mpiexec_cmd.append("-p") + proc = Popen(mpiexec_cmd, stdout=PIPE, stderr=STDOUT, cwd=self.simpath) while True: @@ -348,7 +372,7 @@ def run_parallel(self, exe): # success is when the success message appears # in every process of the parallel simulation if normal_msg in line.lower(): - nr_success = nr_success + 1 + nr_success += 1 if nr_success == self.ncpus: success = True line = line.rstrip("\r\n") @@ -359,7 +383,6 @@ def run_parallel(self, exe): return success, buff - def compare(self): """ Compare the model results @@ -411,7 +434,6 @@ def compare(self): ext = os.path.splitext(file1)[1][1:] if ext.lower() in head_extensions: - # simulation file pth = os.path.join(self.simpath, file1) files1.append(pth) @@ -425,7 +447,6 @@ def compare(self): # Check to see if there is a corresponding compare file if files_cmp is not None: - if file1 + ".cmp" in files_cmp: # compare file idx = files_cmp.index(file1 + ".cmp") diff --git a/autotest/test_par_gwf02.py b/autotest/test_par_gwf02.py index d8681868dae..b94abb29cf5 100644 --- a/autotest/test_par_gwf02.py +++ b/autotest/test_par_gwf02.py @@ -8,7 +8,8 @@ from simulation import TestSimulation # Test for parallel MODFLOW running a simple -# multi-model setup on different partitionings +# multi-model setup with different numbers +# of partitions # # # [M1ny] | ... | ... | [Mnxny] diff --git a/autotest/test_par_gwf_ims_csv.py b/autotest/test_par_gwf_ims_csv.py new file mode 100644 index 00000000000..96408e17756 --- /dev/null +++ b/autotest/test_par_gwf_ims_csv.py @@ -0,0 +1,263 @@ +import os + +import flopy +import numpy as np +import pytest +from framework import TestFramework +from simulation import TestSimulation + +# Test for parallel MODFLOW running on two cpus. +# It contains two coupled models with +# +# 1d: (nlay,nrow,ncol) = (1,1,5), +# +# constant head boundaries left=1.0, right=10.0. +# The result should be a uniform flow field. + +ex = ["par_gwf_csv"] +dis_shape = [(1, 1, 5)] + +# global convenience... +name_left = "leftmodel" +name_right = "rightmodel" + +# solver data +nouter, ninner = 100, 300 +hclose, rclose, relax = 10e-6, 1e-3, 0.97 + + +def get_model(idx, dir): + name = ex[idx] + + # parameters and spd + # tdis + nper = 1 + tdis_rc = [] + for i in range(nper): + tdis_rc.append((1.0, 1, 1)) + + # model spatial discretization + nlay = dis_shape[idx][0] + nrow = dis_shape[idx][1] + ncol = dis_shape[idx][2] + + # cell spacing + delr = 100.0 + delc = 100.0 + area = delr * delc + + # shift + shift_x = 5 * delr + shift_y = 0.0 + + # top/bot of the aquifer + tops = [0.0, -100.0, -200.0, -300.0, -400.0, -500.0] + + # hydraulic conductivity + k11 = 1.0 + + # boundary stress period data + h_left = 1.0 + h_right = 10.0 + + # initial head + h_start = 0.0 + + sim = flopy.mf6.MFSimulation( + sim_name=name, + version="mf6", + exe_name="mf6", + sim_ws=dir, + ) + + tdis = flopy.mf6.ModflowTdis( + sim, time_units="DAYS", nper=nper, perioddata=tdis_rc + ) + + ims = flopy.mf6.ModflowIms( + sim, + print_option="ALL", + csv_outer_output_filerecord=f"{name}.outer.csv", + csv_inner_output_filerecord=f"{name}.inner.csv", + outer_dvclose=hclose, + outer_maximum=nouter, + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=rclose, + linear_acceleration="BICGSTAB", + relaxation_factor=relax, + ) + + # submodel on the left: + left_chd = [ + [(ilay, irow, 0), h_left] + for irow in range(nrow) + for ilay in range(nlay) + ] + chd_spd_left = {0: left_chd} + + gwf = flopy.mf6.ModflowGwf( + sim, + modelname=name_left, + ) + dis = flopy.mf6.ModflowGwfdis( + gwf, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=tops[0], + botm=tops[1 : nlay + 1], + ) + ic = flopy.mf6.ModflowGwfic( + gwf, + strt=h_start, + ) + npf = flopy.mf6.ModflowGwfnpf( + gwf, + icelltype=0, + k=k11, + ) + chd = flopy.mf6.ModflowGwfchd( + gwf, + stress_period_data=chd_spd_left, + ) + oc = flopy.mf6.ModflowGwfoc( + gwf, + head_filerecord=f"{name_left}.hds", + printrecord=[("HEAD", "LAST"), ("BUDGET", "LAST")], + saverecord=[("HEAD", "LAST")], + ) + + # submodel on the right: + right_chd = [ + [(ilay, irow, ncol - 1), h_right] + for irow in range(nrow) + for ilay in range(nlay) + ] + chd_spd_right = {0: right_chd} + + gwf = flopy.mf6.ModflowGwf( + sim, + modelname=name_right, + ) + dis = flopy.mf6.ModflowGwfdis( + gwf, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + xorigin=shift_x, + yorigin=shift_y, + top=tops[0], + botm=tops[1 : nlay + 1], + ) + ic = flopy.mf6.ModflowGwfic( + gwf, + strt=h_start, + ) + npf = flopy.mf6.ModflowGwfnpf( + gwf, + icelltype=0, + k=k11, + ) + chd = flopy.mf6.ModflowGwfchd( + gwf, + stress_period_data=chd_spd_right, + ) + oc = flopy.mf6.ModflowGwfoc( + gwf, + head_filerecord=f"{name_right}.hds", + printrecord=[("HEAD", "LAST"), ("BUDGET", "LAST")], + saverecord=[("HEAD", "LAST")], + ) + + # exchangedata + angldegx = 0.0 + cdist = delr + gwfgwf_data = [ + [ + (ilay, irow, ncol - 1), + (ilay, irow, 0), + 1, + delr / 2.0, + delr / 2.0, + delc, + angldegx, + cdist, + ] + for irow in range(nrow) + for ilay in range(nlay) + ] + gwfgwf = flopy.mf6.ModflowGwfgwf( + sim, + exgtype="GWF6-GWF6", + nexg=len(gwfgwf_data), + exgmnamea=name_left, + exgmnameb=name_right, + exchangedata=gwfgwf_data, + auxiliary=["ANGLDEGX", "CDIST"], + ) + + return sim + + +def build_petsc_db(exdir): + petsc_db_file = os.path.join(exdir, ".petscrc") + with open(petsc_db_file, "w") as petsc_file: + petsc_file.write("-ksp_type cg\n") + petsc_file.write("-pc_type bjacobi\n") + petsc_file.write("-sub_pc_type ilu\n") + petsc_file.write("-sub_pc_factor_levels 2\n") + petsc_file.write(f"-dvclose {hclose}\n") + petsc_file.write(f"-ksp_max_it {nouter}\n") + petsc_file.write("-options_left no\n") + # petsc_file.write("-log_view\n") + + +def build_model(idx, exdir): + sim = get_model(idx, exdir) + build_petsc_db(exdir) + return sim, None + + +def eval_model(sim): + # two coupled models with a uniform flow field, + # here we assert the known head values at the + # cell centers + fpth = os.path.join(sim.simpath, f"{name_left}.hds") + hds = flopy.utils.HeadFile(fpth) + heads_left = hds.get_data().flatten() + fpth = os.path.join(sim.simpath, f"{name_right}.hds") + hds = flopy.utils.HeadFile(fpth) + heads_right = hds.get_data().flatten() + np.testing.assert_array_almost_equal( + heads_left[0:5], [1.0, 2.0, 3.0, 4.0, 5.0] + ) + np.testing.assert_array_almost_equal( + heads_right[0:5], [6.0, 7.0, 8.0, 9.0, 10.0] + ) + + +@pytest.mark.parallel +@pytest.mark.parametrize( + "idx, name", + list(enumerate(ex)), +) +def test_mf6model(idx, name, function_tmpdir, targets): + test = TestFramework() + test.build(build_model, idx, str(function_tmpdir)) + test.run( + TestSimulation( + name=name, + exe_dict=targets, + exfunc=eval_model, + idxsim=idx, + make_comparison=False, + parallel=True, + ncpus=2, + ), + str(function_tmpdir), + ) diff --git a/meson.build b/meson.build index adc2078da00..ed4d6f4e087 100644 --- a/meson.build +++ b/meson.build @@ -42,7 +42,12 @@ if fc_id == 'gcc' '-Wno-maybe-uninitialized', # "Uninitialized" flags produce false positives with allocatables '-Wno-uninitialized', ] - + if not get_option('parallel') + link_args += [ + '-static-libgfortran', + ] + endif + # Options specific to profile if profile == 'release' compile_args += ['-ffpe-summary=overflow', '-ffpe-trap=overflow,zero,invalid'] diff --git a/src/Solution/NumericalSolution.f90 b/src/Solution/NumericalSolution.f90 index c27ca4df333..7618082d69b 100644 --- a/src/Solution/NumericalSolution.f90 +++ b/src/Solution/NumericalSolution.f90 @@ -19,7 +19,7 @@ module NumericalSolutionModule use BaseSolutionModule, only: BaseSolutionType, AddBaseSolutionToList use ListModule, only: ListType use ListsModule, only: basesolutionlist - use InputOutputModule, only: getunit + use InputOutputModule, only: getunit, append_processor_id use NumericalModelModule, only: NumericalModelType, & AddNumericalModelToList, & GetNumericalModelFromList @@ -27,7 +27,8 @@ module NumericalSolutionModule AddNumericalExchangeToList, & GetNumericalExchangeFromList use SparseModule, only: sparsematrix - use SimVariablesModule, only: iout, isim_mode, errmsg + use SimVariablesModule, only: iout, isim_mode, errmsg, & + proc_id, nr_procs use SimStagesModule use BlockParserModule, only: BlockParserType use IMSLinearModule @@ -541,7 +542,7 @@ subroutine sln_ar(this) integer(I4B) :: ifdparam, mxvl, npp integer(I4B) :: ims_lin_type integer(I4B) :: ierr - logical :: isfound, endOfBlock + logical(LGP) :: isfound, endOfBlock integer(I4B) :: ival real(DP) :: rval character(len=*), parameter :: fmtcsvout = & @@ -610,6 +611,9 @@ subroutine sln_ar(this) call this%parser%GetStringCaps(keyword) if (keyword == 'FILEOUT') then call this%parser%GetString(fname) + if (nr_procs > 1) then + call append_processor_id(fname, proc_id) + end if this%icsvouterout = getunit() call openfile(this%icsvouterout, iout, fname, 'CSV_OUTER_OUTPUT', & filstat_opt='REPLACE') @@ -623,6 +627,9 @@ subroutine sln_ar(this) call this%parser%GetStringCaps(keyword) if (keyword == 'FILEOUT') then call this%parser%GetString(fname) + if (nr_procs > 1) then + call append_processor_id(fname, proc_id) + end if this%icsvinnerout = getunit() call openfile(this%icsvinnerout, iout, fname, 'CSV_INNER_OUTPUT', & filstat_opt='REPLACE') @@ -912,17 +919,19 @@ subroutine sln_ar(this) if (ims_lin_type .eq. 1) then this%isymmetric = 1 end if - ! - ! -- incorrect linear solver flag + ! + ! -- petsc linear solver flag else if (this%linmeth == 2) then call this%linear_solver%initialize(this%system_matrix) this%nitermax = this%linear_solver%nitermax this%isymmetric = 0 - ELSE - WRITE (errmsg, '(a)') & + ! + ! -- incorrect linear solver flag + else + write (errmsg, '(a)') & 'Incorrect value for linear solution method specified.' call store_error(errmsg) - END IF + end if ! ! -- write message about matrix symmetry @@ -1374,13 +1383,18 @@ subroutine writeCSVHeader(this) 'solution_inner_dvmax_node' write (this%icsvinnerout, '(*(G0,:,","))', advance='NO') & '', 'solution_inner_drmax', 'solution_inner_drmax_model', & - 'solution_inner_drmax_node', 'solution_inner_alpha' - if (this%imslinear%ilinmeth == 2) then + 'solution_inner_drmax_node' + ! solver items specific to ims solver + if (this%linmeth == 1) then write (this%icsvinnerout, '(*(G0,:,","))', advance='NO') & - '', 'solution_inner_omega' + '', 'solution_inner_alpha' + if (this%imslinear%ilinmeth == 2) then + write (this%icsvinnerout, '(*(G0,:,","))', advance='NO') & + '', 'solution_inner_omega' + end if end if - ! -- check for more than one model - if (this%convnmod > 1) then + ! -- check for more than one model - ims only + if (this%linmeth == 1 .and. this%convnmod > 1) then do im = 1, this%modellist%Count() mp => GetNumericalModelFromList(this%modellist, im) write (this%icsvinnerout, '(*(G0,:,","))', advance='NO') & @@ -2141,12 +2155,14 @@ subroutine csv_convergence_summary(this, iu, totim, kper, kstp, kouter, & call this%sln_get_nodeu(locdr, im, nodeu) write (iu, '(*(G0,:,","))', advance='NO') '', dr, im, nodeu ! - ! -- write acceleration parameters - write (iu, '(*(G0,:,","))', advance='NO') & - '', trim(adjustl(this%caccel(kpos))) + ! -- write ims acceleration parameters + if (this%linmeth == 1) then + write (iu, '(*(G0,:,","))', advance='NO') & + '', trim(adjustl(this%caccel(kpos))) + end if ! - ! -- write information for each model - if (this%convnmod > 1) then + ! -- write information for each model - ims only + if (this%linmeth == 1 .and. this%convnmod > 1) then do j = 1, this%convnmod locdv = this%convlocdv(j, kpos) dv = this%convdvmax(j, kpos) @@ -2378,7 +2394,7 @@ subroutine sln_ls(this, kiter, kstp, kper, in_iter, iptc, ptcf) integer(I4B), intent(inout) :: iptc real(DP), intent(in) :: ptcf ! -- local variables - logical :: lsame + logical(LGP) :: lsame integer(I4B) :: ieq integer(I4B) :: irow integer(I4B) :: itestmat @@ -3143,10 +3159,10 @@ function sln_package_convergence(this, dpak, cpakout, iend) result(ivalue) ! dummy class(NumericalSolutionType) :: this !< NumericalSolutionType instance real(DP), intent(in) :: dpak !< Newton Under-relaxation flag - character(len=LENPAKLOC), intent(in) :: cpakout - integer(I4B), intent(in) :: iend + character(len=LENPAKLOC), intent(in) :: cpakout !< string with package that caused failure + integer(I4B), intent(in) :: iend !< flag indicating if last inner iteration (iend=1) ! local - integer(I4B) :: ivalue !< + integer(I4B) :: ivalue ivalue = 1 if (abs(dpak) > this%dvclose) then ivalue = 0 diff --git a/src/Utilities/InputOutput.f90 b/src/Utilities/InputOutput.f90 index ed27ffd90e6..f81c44b48eb 100644 --- a/src/Utilities/InputOutput.f90 +++ b/src/Utilities/InputOutput.f90 @@ -24,7 +24,8 @@ module InputOutputModule GetFileFromPath, extract_idnum_or_bndname, urdaux, & get_jk, print_format, BuildFixedFormat, & BuildFloatFormat, BuildIntFormat, fseek_stream, & - get_nwords, u9rdcom + get_nwords, u9rdcom, & + append_processor_id contains @@ -256,6 +257,42 @@ subroutine lowcase(word) return end subroutine lowcase + !> @brief Append processor id to a string + !! + !! Subroutine to append the processor id to a string + !! before the file extension (extension is the + !! string after the last '.' in the string. If there + !! is no '.' in the string the processor id is appended + !! to the end of the string. + !! + !< + subroutine append_processor_id(name, proc_id) + ! -- dummy variables + character(len=linelength), intent(inout) :: name !< file name + integer(I4B), intent(in) :: proc_id !< processor id + ! -- local variables + character(len=linelength) :: name_local + character(len=linelength) :: extension_local + integer(I4B) :: ipos0 + integer(I4B) :: ipos1 + ! + name_local = name + call lowcase(name_local) + ipos0 = index(name_local, ".", back=.TRUE.) + ipos1 = len_trim(name) + if (ipos0 > 0) then + write(extension_local, '(a)') name(ipos0:ipos1) + else + ipos0 = ipos1 + extension_local = '' + end if + write(name, '(a,a,i0,a)') & + name(1:ipos0-1), '.p', proc_id, trim(adjustl(extension_local)) + ! + ! -- return + return + end subroutine append_processor_id + !> @brief Create a formatted line !! !! Subroutine to create a formatted line with specified alignment diff --git a/src/mf6core.f90 b/src/mf6core.f90 index f173596e33f..6ec3d061e3e 100644 --- a/src/mf6core.f90 +++ b/src/mf6core.f90 @@ -229,7 +229,7 @@ end subroutine print_info subroutine create_lstfile() use ConstantsModule, only: LINELENGTH use SimVariablesModule, only: proc_id, nr_procs, simlstfile, iout - use InputOutputModule, only: getunit, openfile + use InputOutputModule, only: getunit, openfile, append_processor_id use GenericUtilitiesModule, only: sim_message use VersionModule, only: write_listfile_header character(len=LINELENGTH) :: line @@ -238,7 +238,7 @@ subroutine create_lstfile() iout = getunit() ! if (nr_procs > 1) then - write (simlstfile, '(a,i0,a)') 'mfsim.p', proc_id, '.lst' + call append_processor_id(simlstfile, proc_id) end if ! call openfile(iout, 0, simlstfile, 'LIST', filstat_opt='REPLACE')