From 9a22f639865b9d81716f7fcecdfe33b3e5a352a2 Mon Sep 17 00:00:00 2001 From: Ralf Junghanns Date: Thu, 3 Aug 2023 23:03:21 +0200 Subject: [PATCH] fix(mtlistfile): fix reading MT3D budget (#1899) * fix character span reading tkstp with a regex match * skip particle information reading GW-Mass Budget Co-authored-by: Ralf Junghanns --- autotest/test_listbudget.py | 4 + examples/data/mt3d_test/mt3d_with_adv.list | 330 +++++++++++++++++++++ flopy/utils/mtlistfile.py | 44 ++- 3 files changed, 374 insertions(+), 4 deletions(-) create mode 100644 examples/data/mt3d_test/mt3d_with_adv.list diff --git a/autotest/test_listbudget.py b/autotest/test_listbudget.py index 7f39e8d08..a9d7ce792 100644 --- a/autotest/test_listbudget.py +++ b/autotest/test_listbudget.py @@ -116,6 +116,10 @@ def test_mtlist(example_data_path): mt = MtListBudget(mt_dir / "mcomp.list") df_gw, df_sw = mt.parse(forgive=False, diff=False, start_datetime=None) + mt_dir = example_data_path / "mt3d_test" + mt = MtListBudget(mt_dir / "mt3d_with_adv.list") + df_gw, df_sw = mt.parse(forgive=False, diff=False, start_datetime=None) + mt_dir = example_data_path / "mt3d_test" mt = MtListBudget(mt_dir / "CrnkNic.mt3d.list") df_gw, df_sw = mt.parse(forgive=False, diff=True, start_datetime=None) diff --git a/examples/data/mt3d_test/mt3d_with_adv.list b/examples/data/mt3d_test/mt3d_with_adv.list new file mode 100644 index 000000000..daf39ef40 --- /dev/null +++ b/examples/data/mt3d_test/mt3d_with_adv.list @@ -0,0 +1,330 @@ + LISTING FILE: mt.list + UNIT 16 + + OPENING mt3d_link.ftl + FILE TYPE:FTL UNIT 10 + + OPENING mt.btn + FILE TYPE:BTN UNIT 31 + + OPENING mt.adv + FILE TYPE:ADV UNIT 32 + + OPENING mt.dsp + FILE TYPE:DSP UNIT 33 + + OPENING mt.gcg + FILE TYPE:GCG UNIT 35 + + OPENING mt.ssm + FILE TYPE:SSM UNIT 34 + + OPENING mt.rct + FILE TYPE:RCT UNIT 36 + + +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + + + + MT3DMS + + + A Modular 3D Multi-Species Transport Model + + + For Simulation of Advection, Dispersion and Chemical Reactions + + + of Contaminants in Groundwater Systems + + + + + +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + ----- + | M T | ## BTN for MT3DMS, generated by Flopy. + | 3 D | ## + ----- + THE TRANSPORT MODEL CONSISTS OF 1 LAYER(S) 31 ROW(S) 31 COLUMN(S) + NUMBER OF STRESS PERIOD(S) FOR TRANSPORT SIMULATION = 1 + NUMBER OF ALL COMPONENTS INCLUDED IN SIMULATION = 1 + NUMBER OF MOBILE COMPONENTS INCLUDED IN SIMULATION = 1 + UNIT FOR TIME IS D ; UNIT FOR LENGTH IS M ; UNIT FOR MASS IS KG + OPTIONAL PACKAGES INCLUDED IN CURRENT SIMULATION: + o ADV ON UNIT 32 + o DSP ON UNIT 33 + o SSM ON UNIT 34 + o RCT ON UNIT 36 + o GCG ON UNIT 35 + + BTN5 -- BASIC TRANSPORT PACKAGE, VERSION 5, FEBRUARY 2010, INPUT READ FROM UNIT 31 + 15500 ELEMENTS OF THE X ARRAY USED BY THE BTN PACKAGE + 962 ELEMENTS OF THE IX ARRAY USED BY THE BTN PACKAGE + + FMI5 -- FLOW MODEL INTERFACE PACKAGE, VERSION 5, FEBRUARY 2010, INPUT READ FROM UNIT 10 + FLOW MODEL IS STEADY-STATE + + ADV5 -- ADVECTION PACKAGE, VERSION 5, FEBRUARY 2010, INPUT READ FROM UNIT 32 + ADVECTION IS SOLVED WITH THE HYBRID [MOC]/[MMOC] SCHEME + COURANT NUMBER ALLOWED IN SOLVING THE ADVECTION TERM = 0.750 + MAXIMUM NUMBER OF MOVING PARTICLES ALLOWED = 800000 + 3200000 ELEMENTS OF THE X ARRAY USED BY THE ADV PACKAGE + 1600961 ELEMENTS OF THE IX ARRAY USED BY THE ADV PACKAGE + + DSP5 -- DISPERSION PACKAGE, VERSION 5, FEBRUARY 2010, INPUT READ FROM UNIT 33 + 10573 ELEMENTS OF THE X ARRAY USED BY THE DSP PACKAGE + 0 ELEMENTS OF THE IX ARRAY USED BY THE DSP PACKAGE + + SSM5 -- SINK & SOURCE MIXING PACKAGE, VERSION 5, FEBRUARY 2010, INPUT READ FROM UNIT 34 + HEADER LINE OF THE SSM PACKAGE INPUT FILE: + T F F F F F F F F F F F F F F F + MAJOR STRESS COMPONENTS PRESENT IN THE FLOW MODEL: + o WELL [WEL] + MAXIMUM NUMBER OF POINT SINKS/SOURCES = 121 + 1573 ELEMENTS OF THE X ARRAY USED BY THE SSM PACKAGE + 0 ELEMENTS OF THE IX ARRAY BY THE SSM PACKAGE + + RCT5 -- CHEMICAL REACTION PACKAGE, VERSION 5, FEBRUARY 2010, INPUT READ FROM UNIT 36 + NO SORPTION [OR DUAL-DOMAIN MODEL] IS SIMULATED + NO FIRST-ORDER RATE REACTION IS SIMULATED + REACTION COEFFICIENTS ASSIGNED CELL-BY-CELL + INITIAL SORBED/IMMOBILE PHASE CONCENTRATION ASSIGNED BY DEFAULT + 0 ELEMENTS OF THE X ARRAY USED BY THE RCT PACKAGE + 0 ELEMENTS OF THE IX ARRAY USED BY THE RCT PACKAGE + + GCG5 -- GENERALIZED CONJUGATE GRADIENT SOLVER PACKAGE, VERSION 5, FEBRUARY 2010 INPUT READ FROM UNIT 35 + MAXIMUM OF 1 OUTER ITERATIONS + AND 50 INNER ITERATIONS ALLOWED FOR CLOSURE + THE PRECONDITIONING TYPE SELECTED IS MODIFIED INCOMPLETE CHOLESKY (MIC). + DISPERSION CROSS TERMS LUMPED INTO RIGHT-HAND-SIDE + 21192 ELEMENTS OF THE X ARRAY USED BY THE GCG PACKAGE + 150 ELEMENTS OF THE IX ARRAY USED BY THE GCG PACKAGE + + .......................................... + ELEMENTS OF THE X ARRAY USED = 3248839 + ELEMENTS OF THE IX ARRAY USED = 1602074 + .......................................... + + LAYER NUMBER AQUIFER TYPE + ------------ ------------ + 1 0 + + WIDTH ALONG ROWS (DELR) READ ON UNIT 31 USING FORMAT: " (31E15.6) " + ----------------------------------------------------------------------------- + + WIDTH ALONG COLS (DELC) READ ON UNIT 31 USING FORMAT: " (31E15.6) " + ----------------------------------------------------------------------------- + TOP ELEV. OF 1ST LAYER = 10.00000 + + CELL THICKNESS (DZ) FOR LAYER 1 READ ON UNIT 31 USING FORMAT: " (31E15.6) " + ------------------------------------------------------------------------------------------ + POROSITY = 0.3000000 FOR LAYER 1 + CONCN. BOUNDARY ARRAY = 1 FOR LAYER 1 + INITIAL CONC.: COMP. 01 = 0.000000 FOR LAYER 1 + + VALUE INDICATING INACTIVE CONCENTRATION CELLS = 0.1000000E+31 + MINIMUM SATURATED THICKNESS [THKMIN] ALLOWED = 0.0100 OF TOTAL CELL THICKNESS + + + OUTPUT CONTROL OPTIONS + ---------------------- + + DO NOT PRINT CELL CONCENTRATION + DO NOT PRINT PARTICLE NUMBER IN EACH CELL + DO NOT PRINT RETARDATION FACTOR + DO NOT PRINT DISPERSION COEFFICIENT + SAVE DISSOLVED PHASE CONCENTRATIONS IN UNFORMATTED FILES [MT3Dnnn.UCN] + FOR EACH SPECIES ON UNITS 201 AND ABOVE + + NUMBER OF TIMES AT WHICH SIMULATION RESULTS ARE SAVED = 0 + + NUMBER OF OBSERVATION POINTS = 0 + + SAVE ONE-LINE SUMMARY OF MASS BUDGETS IN FILES [MT3Dnnn.MAS] + FOR EACH SPECIES ON UNITS 601 AND ABOVE, EVERY 1 TRANSPORT STEPS + + MAXIMUM LENGTH ALONG THE X (J) AXIS = 311.0040 + MAXIMUM LENGTH ALONG THE Y (I) AXIS = 313.7939 + MAXIMUM LENGTH ALONG THE Z (K) AXIS = 1.000000 + + + ADVECTION SOLUTION OPTIONS + -------------------------- + + ADVECTION IS SOLVED WITH THE HYBRID [MOC]/[MMOC] SCHEME + COURANT NUMBER ALLOWED IN SOLVING THE ADVECTION TERM = 0.750 + MAXIMUM NUMBER OF MOVING PARTICLES ALLOWED = 800000 + METHOD FOR PARTICLE TRACKING IS [MIXED ORDER] + CONCENTRATION WEIGHTING FACTOR [WD] = 0.500 + THE CONCENTRATION GRADIENT CONSIDERED NEGLIGIBLE [DCEPS] = 0.1000000E-04 + INITIAL PARTICLES ARE PLACED ON 2 VERTICAL PLANE(S) WITHIN CELL BLOCK + PARTICLE NUMBER PER CELL IF DCCELL =< DCEPS = 10 + PARTICLE NUMBER PER CELL IF DCCELL > DCEPS = 40 + MINIMUM PARTICLE NUMBER ALLOWD PER CELL = 5 + MAXIMUM PARTICLE NUMBER ALLOWD PER CELL = 80 + MULTIPLIER OF PARTICLE NUMBER AT SOURCE = 1.00 + SCHEME FOR CONCENTRATION INTERPOLATION IS [LINEAR] + PARTICLES FOR APPROXIMATING A SINK CELL IN THE [MMOC] SCHEME + ARE PLACED RANDOMLY WITHIN CELL BLOCK + NUMBER OF PARTICLES USED TO APPROXIMATE A SINK CELL IN THE [MMOC] SCHEME = 15 + CRITICAL CONCENTRATION GRADIENT USED IN THE "HMOC" SCHEME [DCHMOC] = 0.1000E-03 + THE "MOC" SOLUTION IS USED WHEN DCCELL > DCHMOC + THE "MMOC" SOLUTION IS USED WHEN DCCELL =< DCHMOC + + + DISPERSION INPUT PARAMETERS + --------------------------- + + LONG. DISPERSIVITY (AL) = 0.1000000E-01 FOR LAYER 1 + H. TRANS./LONG. DISP. = 0.1000000 + V. TRANS./LONG. DISP. = 0.1000000E-01 + DIFFUSION COEFFICIENT = 0.1000000E-08 + + + SORPTION AND 1ST/0TH ORDER REACTION PARAMETERS + ---------------------------------------------- + + + + + SOLUTION BY THE GENERALIZED CONJUGATE GRADIENT METHOD + ----------------------------------------------------- + MAXIMUM OUTER ITERATIONS ALLOWED FOR CLOSURE = 1 + MAXIMUM INNER ITERATIONS ALLOWED FOR CLOSURE = 50 + PRECONDITIONING TYPE SELECTED = 3 + ACCELERATION PARAMETER = 1.0000 + CONCENTRATION CHANGE CRITERION FOR CLOSURE = 0.10000E-04 + GCG CONCENTRATION CHANGE PRINTOUT INTERVAL = 999 + + + ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + STRESS PERIOD NO. 001 + ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + + LENGTH OF CURRENT STRESS PERIOD = 26.00000 + NUMBER OF TIME STEPS FOR CURRENT STRESS PERIOD = 1 + TIME STEP MULTIPLIER USED IN FLOW SOLUTION = 1.000000 + + ***Type of Transport Simulation is TRANSIENT + + USER-SPECIFIED TRANSPORT STEPSIZE = 0.000000 D + MAXIMUM NUMBER OF TRANSPORT STEPS ALLOWED IN ONE FLOW TIME STEP = 50000 + MULTIPLIER FOR SUCCESSIVE TRANSPORT STEPS [USED IN IMPLICIT SCHEMES] = 1.000 + MAXIMUM TRANSPORT STEP SIZE [USED IN IMPLICIT SCHEMES] = 0.000000 D + + NO LAYER ROW COLUMN CONCENTRATION TYPE COMPONENT + 1 1 16 16 50.00000 WELL 1 + + + ================================================ + TIME STEP NO. 001 + ================================================ + + FROM TIME = 0.0000 TO 26.000 + + + "THKSAT " FLOW TERMS FOR TIME STEP 1, STRESS PERIOD 1 READ UNFORMATTED ON UNIT 10 + -------------------------------------------------------------------------------------------- + + "QXX " FLOW TERMS FOR TIME STEP 1, STRESS PERIOD 1 READ UNFORMATTED ON UNIT 10 + -------------------------------------------------------------------------------------------- + + "QYY " FLOW TERMS FOR TIME STEP 1, STRESS PERIOD 1 READ UNFORMATTED ON UNIT 10 + -------------------------------------------------------------------------------------------- + + MAXIMUM STEPSIZE DURING WHICH ANY PARTICLE CANNOT MOVE MORE THAN ONE CELL + = 1.666 (WHEN MIN. R.F.=1) AT K= 1, I= 16, J= 15 + + MAXIMUM STEPSIZE WHICH MEETS STABILITY CRITERION OF THE ADVECTION TERM + (FOR PURE FINITE-DIFFERENCE OPTION, MIXELM=0) + = 0.6094 (WHEN MIN. R.F.=1) AT K= 1, I= 16, J= 16 + + "CNH " FLOW TERMS FOR TIME STEP 1, STRESS PERIOD 1 READ UNFORMATTED ON UNIT 10 + -------------------------------------------------------------------------------------------- + + "WEL " FLOW TERMS FOR TIME STEP 1, STRESS PERIOD 1 READ UNFORMATTED ON UNIT 10 + -------------------------------------------------------------------------------------------- + + + TOTAL NUMBER OF POINT SOURCES/SINKS PRESENT IN THE FLOW MODEL = 121 + + MAXIMUM STEPSIZE WHICH MEETS STABILITY CRITERION OF THE SINK & SOURCE TERM + = 0.3047 (WHEN MIN. R.F.=1) AT K= 1, I= 16, J= 16 + + MAXIMUM STEPSIZE WHICH MEETS STABILITY CRITERION OF THE DISPERSION TERM + = 307.1 (WHEN MIN. R.F.=1) AT K= 1, I= 16, J= 16 + + + 1 CALLS TO GCG PACKAGE FOR TRANSPORT TIME STEP 1 IN FLOW TIME STEP 1 STRESS PERIOD 1 + 2 TOTAL ITERATIONS + 1 CALLS TO GCG PACKAGE FOR TRANSPORT TIME STEP 2 IN FLOW TIME STEP 1 STRESS PERIOD 1 + 2 TOTAL ITERATIONS + 1 CALLS TO GCG PACKAGE FOR TRANSPORT TIME STEP 3 IN FLOW TIME STEP 1 STRESS PERIOD 1 + 2 TOTAL ITERATIONS + 1 CALLS TO GCG PACKAGE FOR TRANSPORT TIME STEP 4 IN FLOW TIME STEP 1 STRESS PERIOD 1 + 2 TOTAL ITERATIONS + 1 CALLS TO GCG PACKAGE FOR TRANSPORT TIME STEP 5 IN FLOW TIME STEP 1 STRESS PERIOD 1 + 2 TOTAL ITERATIONS + 1 CALLS TO GCG PACKAGE FOR TRANSPORT TIME STEP 6 IN FLOW TIME STEP 1 STRESS PERIOD 1 + 2 TOTAL ITERATIONS + 1 CALLS TO GCG PACKAGE FOR TRANSPORT TIME STEP 7 IN FLOW TIME STEP 1 STRESS PERIOD 1 + 2 TOTAL ITERATIONS + 1 CALLS TO GCG PACKAGE FOR TRANSPORT TIME STEP 8 IN FLOW TIME STEP 1 STRESS PERIOD 1 + 2 TOTAL ITERATIONS + 1 CALLS TO GCG PACKAGE FOR TRANSPORT TIME STEP 9 IN FLOW TIME STEP 1 STRESS PERIOD 1 + 2 TOTAL ITERATIONS + 1 CALLS TO GCG PACKAGE FOR TRANSPORT TIME STEP 10 IN FLOW TIME STEP 1 STRESS PERIOD 1 + 2 TOTAL ITERATIONS + 1 CALLS TO GCG PACKAGE FOR TRANSPORT TIME STEP 11 IN FLOW TIME STEP 1 STRESS PERIOD 1 + 2 TOTAL ITERATIONS + 1 CALLS TO GCG PACKAGE FOR TRANSPORT TIME STEP 12 IN FLOW TIME STEP 1 STRESS PERIOD 1 + 2 TOTAL ITERATIONS + 1 CALLS TO GCG PACKAGE FOR TRANSPORT TIME STEP 13 IN FLOW TIME STEP 1 STRESS PERIOD 1 + 2 TOTAL ITERATIONS + 1 CALLS TO GCG PACKAGE FOR TRANSPORT TIME STEP 14 IN FLOW TIME STEP 1 STRESS PERIOD 1 + 2 TOTAL ITERATIONS + 1 CALLS TO GCG PACKAGE FOR TRANSPORT TIME STEP 15 IN FLOW TIME STEP 1 STRESS PERIOD 1 + 2 TOTAL ITERATIONS + 1 CALLS TO GCG PACKAGE FOR TRANSPORT TIME STEP 16 IN FLOW TIME STEP 1 STRESS PERIOD 1 + 2 TOTAL ITERATIONS + 1 CALLS TO GCG PACKAGE FOR TRANSPORT TIME STEP 17 IN FLOW TIME STEP 1 STRESS PERIOD 1 + 2 TOTAL ITERATIONS + 1 CALLS TO GCG PACKAGE FOR TRANSPORT TIME STEP 18 IN FLOW TIME STEP 1 STRESS PERIOD 1 + 2 TOTAL ITERATIONS + 1 CALLS TO GCG PACKAGE FOR TRANSPORT TIME STEP 19 IN FLOW TIME STEP 1 STRESS PERIOD 1 + 2 TOTAL ITERATIONS + 1 CALLS TO GCG PACKAGE FOR TRANSPORT TIME STEP 20 IN FLOW TIME STEP 1 STRESS PERIOD 1 + 2 TOTAL ITERATIONS + 1 CALLS TO GCG PACKAGE FOR TRANSPORT TIME STEP 21 IN FLOW TIME STEP 1 STRESS PERIOD 1 + 2 TOTAL ITERATIONS + MAXIMUM CONCENTRATION CHANGES FOR EACH ITERATION: + MAX. CHANGE LAYER,ROW,COL MAX. CHANGE LAYER,ROW,COL MAX. CHANGE LAYER,ROW,COL MAX. CHANGE LAYER,ROW,COL MAX. CHANGE LAYER,ROW,COL + ------------------------------------------------------------------------------------------------------------------------------------ + 0.1276E-03 ( 1, 14, 12) 0.7523E-36 ( 1, 3, 30) + + + + >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>FOR COMPONENT NO. 01<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + ------------------------------------------- + TRANSPORT STEP NO. 21 + ------------------------------------------- + + TOTAL ELAPSED TIME SINCE BEGINNING OF SIMULATION = 26.00000 D + ..................................................................... + + TOTAL PARTICLES USED IN THE CURRENT STEP = 4264 + PARTICLES ADDED AT BEGINNING OF THE STEP = 256 + PARTICLES REMOVED AT END OF LAST STEP = 0 + + CUMMULATIVE MASS BUDGETS AT END OF TRANSPORT STEP 21, TIME STEP 1, STRESS PERIOD 1 + ------------------------------------------------------------------------------------------ + + IN OUT + ---------------- ---------------- + CONSTANT CONCENTRATION: 0.000000 0.000000 + CONSTANT HEAD: 0.000000 -0.1140626E-09 + WELLS: 130000.0 0.000000 + 1ST/0TH ORDER REACTION: 0.000000 0.000000 + MASS STORAGE (SOLUTE): 43.12313 -122520.2 + --------------------------------------------------------------------------- + [TOTAL]: 130043.1 KG -122520.2 KG + + NET (IN - OUT): 7522.922 + DISCREPANCY (PERCENT): 5.957257 + ----- + | M T | + | 3 D | END OF MODEL OUTPUT + ----- diff --git a/flopy/utils/mtlistfile.py b/flopy/utils/mtlistfile.py index 9f70c72b4..af682bc7a 100644 --- a/flopy/utils/mtlistfile.py +++ b/flopy/utils/mtlistfile.py @@ -3,13 +3,12 @@ mt3d(usgs) run. Also includes support for SFT budget. """ +import re import warnings import numpy as np import pandas as pd -from ..utils import import_optional_dependency - class MtListBudget: """ @@ -24,7 +23,6 @@ class MtListBudget: Examples -------- >>> mt_list = MtListBudget("my_mt3d.list") - >>> incremental, cumulative = mt_list.get_budget() >>> gw_df, sw_df = mt_list.parse(start_datetime="10-21-2015") """ @@ -48,6 +46,7 @@ def __init__(self, file_name): self.time_key = line.lower() line = "TRANSPORT TIME STEP" self.tkstp_key = line.lower() + self.particles_key = "TOTAL PARTICLES USED IN THE CURRENT STEP" return @@ -110,7 +109,18 @@ def parse( else: self._parse_sw(f, line) elif self.tkstp_key in line: - self.tkstp_overflow = int(line[51:58]) + try: + self.tkstp_overflow = ( + self._extract_number_between_strings( + line, self.tkstp_key, "in" + ) + ) + except Exception as e: + warnings.warn( + "error parsing TKSTP key " + f"starting on line {self.lcount}: {e!s}" + ) + break if len(self.gw_data) == 0: raise Exception("no groundwater budget info found...") @@ -257,6 +267,15 @@ def _parse_gw(self, f, line): line = self._readline(f) if line is None: raise Exception("EOF while reading from totim to time step") + + if self.particles_key.lower() in line: + for _ in range(4): + line = self._readline(f) + if line is None: + raise Exception( + "EOF while reading from time step to particles" + ) + try: kper = int(line[-6:-1]) kstp = int(line[-26:-21]) @@ -472,3 +491,20 @@ def _add_to_sw_data(self, inout, item, cval, fval, comp): if iitem not in self.sw_data.keys(): self.sw_data[iitem] = [] self.sw_data[iitem].append(val) + + @staticmethod + def _extract_number_between_strings( + input_string, start_string, end_string + ): + pattern = ( + rf"{re.escape(start_string)}\s*(\d+)\s*{re.escape(end_string)}" + ) + match = re.search(pattern, input_string) + + if match: + extracted_number = int(match.group(1)) + return extracted_number + else: + raise Exception( + f"Error extracting number between {start_string} and {end_string} in {input_string}" + )