diff --git a/tests/test_run.py b/tests/test_run.py index 8f9586c1..186308f5 100644 --- a/tests/test_run.py +++ b/tests/test_run.py @@ -16,6 +16,7 @@ ], ) def test_run(exec, workers, docker, wsl, cmd): + toughio._run._check_exec = False status = toughio.run( exec, {}, diff --git a/toughio/VERSION b/toughio/VERSION index f0df1f7d..cd99d386 100644 --- a/toughio/VERSION +++ b/toughio/VERSION @@ -1 +1 @@ -1.13.2 \ No newline at end of file +1.14.0 \ No newline at end of file diff --git a/toughio/_io/input/_helpers.py b/toughio/_io/input/_helpers.py index 9b4073f4..78230a4b 100644 --- a/toughio/_io/input/_helpers.py +++ b/toughio/_io/input/_helpers.py @@ -61,6 +61,8 @@ def read(filename, file_format=None, **kwargs): Other Parameters ---------------- + blocks : list of str or None, optional, default None + Only if ``file_format = "tough"``. Blocks to read. If None, all blocks are read. label_length : int or None, optional, default None Only if ``file_format = "tough"``. Number of characters in cell labels. n_variables : int or None, optional, default None diff --git a/toughio/_io/input/tough/_read.py b/toughio/_io/input/tough/_read.py index f09d88d5..10362be5 100644 --- a/toughio/_io/input/tough/_read.py +++ b/toughio/_io/input/tough/_read.py @@ -17,7 +17,14 @@ } -def read(filename, label_length=None, n_variables=None, eos=None, simulator="tough"): +def read( + filename, + blocks=None, + label_length=None, + n_variables=None, + eos=None, + simulator="tough", +): """ Read TOUGH input file. @@ -25,6 +32,8 @@ def read(filename, label_length=None, n_variables=None, eos=None, simulator="tou ---------- filename : str, pathlike or buffer Input file name or buffer. + blocks : list of str or None, optional, default None + Blocks to read. If None, all blocks are read. label_length : int or None, optional, default None Number of characters in cell labels. n_variables : int or None, optional, default None @@ -46,39 +55,46 @@ def read(filename, label_length=None, n_variables=None, eos=None, simulator="tou raise ValueError() with open_file(filename, "r") as f: - out = read_buffer(f, label_length, n_variables, eos, simulator) + out = read_buffer(f, blocks, label_length, n_variables, eos, simulator) return out -def read_buffer(f, label_length, n_variables, eos, simulator="tough"): +def read_buffer(f, block_stack, label_length, n_variables, eos, simulator="tough"): """Read TOUGH input file.""" from ._common import blocks + # Block filters + block_stack = block_stack if block_stack is not None else blocks + block_stack = set(block_stack) + parameters = {} flag = False # Title - title = [] - while True: - if len(title) >= 100: - raise ValueError() + if "TITLE" in block_stack: + block_stack.remove("TITLE") - line = f.readline().strip() + title = [] + while True: + if len(title) >= 100: + raise ValueError() - if line[:5].rstrip().upper() not in blocks: - title.append(line) + line = f.readline().strip() - else: - break + if line[:5].rstrip().upper() not in blocks: + title.append(line) + + else: + break - if title: - title = title[0] if len(title) == 1 else title + if title: + title = title[0] if len(title) == 1 else title - if title: - parameters["title"] = title + if title: + parameters["title"] = title - f.seek(0) + f.seek(0) # Loop over blocks # Some blocks (INCON, INDOM, PARAM) need to rewind to previous line but tell and seek are disabled by next @@ -87,13 +103,16 @@ def read_buffer(f, label_length, n_variables, eos, simulator="tough"): try: for line in fiter: - if line.startswith("DIMEN"): + if line.startswith("DIMEN") and "DIMEN" in block_stack: + block_stack.remove("DIMEN") parameters.update(_read_dimen(fiter)) - elif line.startswith("ROCKS"): + elif line.startswith("ROCKS") and "ROCKS" in block_stack: + block_stack.remove("ROCKS") parameters.update(_read_rocks(fiter, simulator)) - elif line.startswith("RPCAP"): + elif line.startswith("RPCAP") and "RPCAP" in block_stack: + block_stack.remove("RPCAP") rpcap = _read_rpcap(fiter) if "default" in parameters: @@ -102,43 +121,53 @@ def read_buffer(f, label_length, n_variables, eos, simulator="tough"): else: parameters["default"] = rpcap - elif line.startswith("REACT"): + elif line.startswith("REACT") and "REACT" in block_stack: + block_stack.remove("REACT") react = _read_react(fiter) + if "react" in parameters: parameters["react"].update(react["react"]) else: parameters.update(react) - elif line.startswith("FLAC"): + elif line.startswith("FLAC") and "FLAC" in block_stack: + block_stack.remove("FLAC") flac = _read_flac(fiter, parameters["rocks"]) parameters["flac"] = flac["flac"] for k, v in flac["rocks"].items(): parameters["rocks"][k].update(v) - elif line.startswith("CHEMP"): + elif line.startswith("CHEMP") and "CHEMP" in block_stack: + block_stack.remove("CHEMP") parameters.update(_read_chemp(fiter)) - elif line.startswith("NCGAS"): + elif line.startswith("NCGAS") and "NCGAS" in block_stack: + block_stack.remove("NCGAS") parameters.update(_read_ncgas(fiter)) - elif line.startswith("MULTI"): + elif line.startswith("MULTI") and "MULTI" in block_stack: + block_stack.remove("MULTI") parameters.update(_read_multi(fiter)) if not n_variables: n_variables = parameters["n_component"] + 1 - elif line.startswith("SOLVR"): + elif line.startswith("SOLVR") and "SOLVR" in block_stack: + block_stack.remove("SOLVR") parameters.update(_read_solvr(fiter)) - elif line.startswith("INDEX"): + elif line.startswith("INDEX") and "INDEX" in block_stack: + block_stack.remove("INDEX") parameters["index"] = True - elif line.startswith("START"): + elif line.startswith("START") and "START" in block_stack: + block_stack.remove("START") parameters["start"] = True - elif line.startswith("PARAM"): + elif line.startswith("PARAM") and "PARAM" in block_stack: + block_stack.remove("PARAM") param, n_variables = _read_param(fiter, n_variables, eos) parameters["options"] = param["options"] parameters["extra_options"] = param["extra_options"] @@ -149,68 +178,86 @@ def read_buffer(f, label_length, n_variables, eos, simulator="tough"): else: parameters["default"] = param["default"] - elif line.startswith("SELEC"): + elif line.startswith("SELEC") and "SELEC" in block_stack: + block_stack.remove("SELEC") parameters.update(_read_selec(fiter)) - elif line.startswith("INDOM"): + elif line.startswith("INDOM") and "INDOM" in block_stack: + block_stack.remove("INDOM") indom, n_variables = _read_indom(fiter, n_variables, eos) for k, v in indom["rocks"].items(): parameters["rocks"][k].update(v) - elif line.startswith("MOMOP"): + elif line.startswith("MOMOP") and "MOMOP" in block_stack: + block_stack.remove("MOMOP") parameters.update(_read_momop(fiter)) - elif line.startswith("TIMES"): + elif line.startswith("TIMES") and "TIMES" in block_stack: + block_stack.remove("TIMES") parameters.update(_read_times(fiter)) - elif line.startswith("HYSTE"): + elif line.startswith("HYSTE") and "HYSTE" in block_stack: + block_stack.remove("HYSTE") parameters.update(_read_hyste(fiter)) - elif line.startswith("FOFT"): + elif line.startswith("FOFT") and "FOFT" in block_stack: + block_stack.remove("FOFT") oft, label_length = _read_oft(fiter, "FOFT", label_length) parameters.update(oft) - elif line.startswith("COFT"): + elif line.startswith("COFT") and "COFT" in block_stack: + block_stack.remove("COFT") oft, label_length = _read_oft(fiter, "COFT", label_length) parameters.update(oft) - elif line.startswith("GOFT"): + elif line.startswith("GOFT") and "GOFT" in block_stack: + block_stack.remove("GOFT") oft, label_length = _read_oft(fiter, "GOFT", label_length) parameters.update(oft) - elif line.startswith("ROFT"): + elif line.startswith("ROFT") and "ROFT" in block_stack: + block_stack.remove("ROFT") parameters.update(_read_roft(fiter)) - elif line.startswith("GENER"): + elif line.startswith("GENER") and "GENER" in block_stack: + block_stack.remove("GENER") gener, flag, label_length = _read_gener(fiter, label_length, simulator) parameters.update(gener) if flag: break - elif line.startswith("TIMBC"): + elif line.startswith("TIMBC") and "TIMBC" in block_stack: + block_stack.remove("TIMBC") parameters.update(_read_timbc(fiter)) - elif line.startswith("DIFFU"): + elif line.startswith("DIFFU") and "DIFFU" in block_stack: + block_stack.remove("DIFFU") parameters.update(_read_diffu(fiter)) - elif line.startswith("OUTPT"): + elif line.startswith("OUTPT") and "OUTPT" in block_stack: + block_stack.remove("OUTPT") outpt = _read_outpt(fiter) + if "react" in parameters: parameters["react"].update(outpt["react"]) + else: parameters.update(outpt) - elif line.startswith("OUTPU"): + elif line.startswith("OUTPU") and "OUTPU" in block_stack: + block_stack.remove("OUTPU") parameters.update(_read_outpu(fiter)) - elif line.startswith("ELEME"): + elif line.startswith("ELEME") and "ELEME" in block_stack: + block_stack.remove("ELEME") eleme, label_length = _read_eleme(fiter, label_length) parameters.update(eleme) parameters["coordinates"] = False - elif line.startswith("COORD"): + elif line.startswith("COORD") and "COORD" in block_stack: + block_stack.remove("COORD") coord = _read_coord(fiter) for k, v in zip(parameters["elements"], coord): @@ -218,14 +265,16 @@ def read_buffer(f, label_length, n_variables, eos, simulator="tough"): parameters["coordinates"] = True - elif line.startswith("CONNE"): + elif line.startswith("CONNE") and "CONNE" in block_stack: + block_stack.remove("CONNE") conne, flag, label_length = _read_conne(fiter, label_length) parameters.update(conne) if flag: break - elif line.startswith("INCON"): + elif line.startswith("INCON") and "INCON" in block_stack: + block_stack.remove("INCON") incon, flag, label_length, n_variables = _read_incon( fiter, label_length, n_variables, eos, simulator ) @@ -234,25 +283,35 @@ def read_buffer(f, label_length, n_variables, eos, simulator="tough"): if flag: break - elif line.startswith("MESHM"): + elif line.startswith("MESHM") and "MESHM" in block_stack: + block_stack.remove("MESHM") parameters.update(_read_meshm(fiter)) - elif line.startswith("POISE"): + elif line.startswith("POISE") and "POISE" in block_stack: + block_stack.remove("POISE") poise = _read_poise(fiter) + if "react" in parameters: parameters["react"].update(poise["react"]) + else: parameters.update(poise) - elif line.startswith("NOVER"): + elif line.startswith("NOVER") and "NOVER" in block_stack: + block_stack.remove("NOVER") parameters["nover"] = True - elif line.startswith("ENDCY"): + elif line.startswith("ENDCY") and "ENDCY" in block_stack: + block_stack.remove("ENDCY") end_comments = read_end_comments(fiter) if end_comments: parameters["end_comments"] = end_comments + # Stop reading if block stack is empty + if not block_stack: + break + except: raise ReadError(f"failed to parse line {fiter.count}.") diff --git a/toughio/_mesh/_mesh.py b/toughio/_mesh/_mesh.py index 1b557d10..e755d4ee 100644 --- a/toughio/_mesh/_mesh.py +++ b/toughio/_mesh/_mesh.py @@ -94,6 +94,70 @@ def __repr__(self): return "\n".join(lines) + def __getitem__(self, islice): + """Slice mesh.""" + if np.ndim(islice) == 0: + islice = [islice] + + if np.ndim(islice) != 1: + raise ValueError() + + if isinstance(islice[0], (bool, np.bool_)): + if len(islice) != self.n_cells: + raise ValueError() + + elif isinstance(islice[0], (int, np.int8, np.int16, np.int32, np.int64)): + if np.max(islice) >= self.n_cells: + raise ValueError() + + tmp = islice + islice = np.zeros(self.n_cells, dtype=bool) + islice[tmp] = True + + else: + raise TypeError() + + count = 0 + cell_idx = [] + nodes_to_keep = set() + + for _, cell_data in self.cells: + cell_idx.append([]) + + for j, cell in enumerate(cell_data): + if islice[count]: + cell_idx[-1].append(j) + + for node in cell: + nodes_to_keep.add(node) + + count += 1 + + # Prune orphaned nodes + node_map = {} + point_idx = [] + count = 0 + + for i in range(self.n_points): + if i in nodes_to_keep: + node_map[i] = count + point_idx.append(i) + count += 1 + + return Mesh( + points=self.points[point_idx], + cells=[ + ( + cell_type, + np.array([[node_map[i] for i in cell] for cell in cell_data[idx]]), + ) + for idx, (cell_type, cell_data) in zip(cell_idx, self.cells) + ], + point_data={k: v[point_idx] for k, v in self.point_data.items()}, + cell_data={k: v[islice] for k, v in self.cell_data.items()}, + field_data={k: v for k, v in self.field_data.items()}, + ) + def extrude_to_3d(self, height=1.0, axis=2, inplace=True): """ Convert a 2D mesh to 3D by extruding cells along given axis. diff --git a/toughio/_run.py b/toughio/_run.py index cc90fe0e..4b079b22 100644 --- a/toughio/_run.py +++ b/toughio/_run.py @@ -6,6 +6,8 @@ import subprocess import tempfile +_check_exec = True # Bool to be monkeypatched in tests + def run( exec, @@ -19,6 +21,8 @@ def run( use_temp=False, ignore_patterns=None, silent=False, + petsc_args=None, + docker_args=None, **kwargs, ): """ @@ -48,6 +52,10 @@ def run( If provided, output files that match the glob-style patterns will be discarded. silent : bool, optional, default False If `True`, nothing will be printed to standard output. + petsc_args : list or None, optional, default None + List of arguments passed to PETSc solver (written to `.petscrc`). + docker_args : list or None, optional, default None + List of arguments passed to `docker run` command. Other Parameters ---------------- @@ -102,6 +110,50 @@ def run( exec = str(exec) exec = f"{os.path.expanduser('~')}/{exec[2:]}" if exec.startswith("~/") else exec + if _check_exec: + if not docker: + if shutil.which(exec) is None: + raise RuntimeError(f"executable '{exec}' not found.") + + else: + # Check if Docker is in the PATH + if shutil.which("docker") is None: + raise RuntimeError("Docker executable not found.") + + # Check if Docker daemon is running + status = subprocess.run( + ["docker", "version"], + stdout=subprocess.PIPE, + universal_newlines=True, + ) + + if "Server" not in status.stdout: + raise RuntimeError( + "cannot connect to the Docker daemon. Is the Docker daemon running?" + ) + + # Check if Docker image exists + status = subprocess.run( + ["docker", "images", "-q", docker], + stdout=subprocess.PIPE, + universal_newlines=True, + ) + + if not status.stdout: + raise RuntimeError(f"image '{docker}' not found.") + + # Check if the executable exists inside the image + status = subprocess.run( + ["docker", "run", "--rm", docker, "which", exec], + stdout=subprocess.PIPE, + universal_newlines=True, + ) + + if not status.stdout: + raise RuntimeError( + f"executable '{exec}' not found in Docker image '{docker}'." + ) + # Working directory working_dir = os.getcwd() if working_dir is None else working_dir working_dir = pathlib.Path(working_dir) @@ -141,6 +193,18 @@ def run( ): shutil.copy(filename, simulation_dir / new_filename.name) + # PETSc arguments + petsc_args = petsc_args if petsc_args else [] + + if petsc_args: + with open(simulation_dir / ".petscrc", "w") as f: + for arg in petsc_args: + if arg.startswith("-"): + f.write(f"{arg} ") + + else: + f.write(f"{arg}\n") + # Output filename output_filename = f"{input_filename.stem}.out" @@ -167,22 +231,60 @@ def run( except AttributeError: uid = "" - cmd = f"docker run --rm {uid} -v {cwd}:/shared -w /shared {docker} {cmd}" + docker_args = docker_args if docker_args else [] + docker_args += [ + "--rm", + uid, + "-v", + f"{cwd}:/shared", + "-w", + "/shared", + ] + cmd = f"docker run {' '.join(docker_args)} {docker} {cmd}" # Use WSL if wsl and is_windows: cmd = f"bash -c '{cmd}'" - kwargs = {} - if silent: - kwargs["stdout"] = subprocess.DEVNULL - kwargs["stderr"] = subprocess.STDOUT + # Run simulation + if not silent: + # See + p = subprocess.Popen( + cmd, + shell=True, + cwd=str(simulation_dir), + stderr=subprocess.STDOUT, + stdout=subprocess.PIPE, + universal_newlines=False, + ) - else: - kwargs["stderr"] = subprocess.PIPE - kwargs["universal_newlines"] = True + stdout = [] + cr = False + for line in open(os.dup(p.stdout.fileno()), newline=""): + # Handle carriage return + # newline in open is converting \r\n as \r moving \r at the end of the previous string + # This is only an issue for Spyder + line = f"\r{line}" if cr else line + cr = line.endswith("\r") + + line = line[:-2] if cr else line + print(line, end="", flush=True) + stdout.append(line) + + status = subprocess.CompletedProcess( + args=p.args, + returncode=0, + stdout="".join(stdout), + ) - status = subprocess.run(cmd, shell=True, cwd=str(simulation_dir), **kwargs) + else: + status = subprocess.run( + cmd, + shell=True, + cwd=str(simulation_dir), + stdout=subprocess.DEVNULL, + stderr=subprocess.STDOUT, + ) # Copy files from temporary directory and delete it if use_temp: