diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..09be0ef --- /dev/null +++ b/.gitignore @@ -0,0 +1,207 @@ +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[codz] +*$py.class + +# C extensions +*.so + +# Distribution / packaging +.Python +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +wheels/ +share/python-wheels/ +*.egg-info/ +.installed.cfg +*.egg +MANIFEST + +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ +.tox/ +.nox/ +.coverage +.coverage.* +.cache +nosetests.xml +coverage.xml +*.cover +*.py.cover +.hypothesis/ +.pytest_cache/ +cover/ + +# Translations +*.mo +*.pot + +# Django stuff: +*.log +local_settings.py +db.sqlite3 +db.sqlite3-journal + +# Flask stuff: +instance/ +.webassets-cache + +# Scrapy stuff: +.scrapy + +# Sphinx documentation +docs/_build/ + +# PyBuilder +.pybuilder/ +target/ + +# Jupyter Notebook +.ipynb_checkpoints + +# IPython +profile_default/ +ipython_config.py + +# pyenv +# For a library or package, you might want to ignore these files since the code is +# intended to run in multiple environments; otherwise, check them in: +# .python-version + +# pipenv +# According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control. +# However, in case of collaboration, if having platform-specific dependencies or dependencies +# having no cross-platform support, pipenv may install dependencies that don't work, or not +# install all needed dependencies. +#Pipfile.lock + +# UV +# Similar to Pipfile.lock, it is generally recommended to include uv.lock in version control. +# This is especially recommended for binary packages to ensure reproducibility, and is more +# commonly ignored for libraries. +#uv.lock + +# poetry +# Similar to Pipfile.lock, it is generally recommended to include poetry.lock in version control. +# This is especially recommended for binary packages to ensure reproducibility, and is more +# commonly ignored for libraries. +# https://python-poetry.org/docs/basic-usage/#commit-your-poetrylock-file-to-version-control +#poetry.lock +#poetry.toml + +# pdm +# Similar to Pipfile.lock, it is generally recommended to include pdm.lock in version control. +# pdm recommends including project-wide configuration in pdm.toml, but excluding .pdm-python. +# https://pdm-project.org/en/latest/usage/project/#working-with-version-control +#pdm.lock +#pdm.toml +.pdm-python +.pdm-build/ + +# pixi +# Similar to Pipfile.lock, it is generally recommended to include pixi.lock in version control. +#pixi.lock +# Pixi creates a virtual environment in the .pixi directory, just like venv module creates one +# in the .venv directory. It is recommended not to include this directory in version control. +.pixi + +# PEP 582; used by e.g. github.com/David-OConnor/pyflow and github.com/pdm-project/pdm +__pypackages__/ + +# Celery stuff +celerybeat-schedule +celerybeat.pid + +# SageMath parsed files +*.sage.py + +# Environments +.env +.envrc +.venv +env/ +venv/ +ENV/ +env.bak/ +venv.bak/ + +# Spyder project settings +.spyderproject +.spyproject + +# Rope project settings +.ropeproject + +# mkdocs documentation +/site + +# mypy +.mypy_cache/ +.dmypy.json +dmypy.json + +# Pyre type checker +.pyre/ + +# pytype static type analyzer +.pytype/ + +# Cython debug symbols +cython_debug/ + +# PyCharm +# JetBrains specific template is maintained in a separate JetBrains.gitignore that can +# be found at https://github.com/github/gitignore/blob/main/Global/JetBrains.gitignore +# and can be added to the global gitignore or merged into this file. For a more nuclear +# option (not recommended) you can uncomment the following to ignore the entire idea folder. +#.idea/ + +# Abstra +# Abstra is an AI-powered process automation framework. +# Ignore directories containing user credentials, local state, and settings. +# Learn more at https://abstra.io/docs +.abstra/ + +# Visual Studio Code +# Visual Studio Code specific template is maintained in a separate VisualStudioCode.gitignore +# that can be found at https://github.com/github/gitignore/blob/main/Global/VisualStudioCode.gitignore +# and can be added to the global gitignore or merged into this file. However, if you prefer, +# you could uncomment the following to ignore the entire vscode folder +# .vscode/ + +# Ruff stuff: +.ruff_cache/ + +# PyPI configuration file +.pypirc + +# Cursor +# Cursor is an AI-powered code editor. `.cursorignore` specifies files/directories to +# exclude from AI features like autocomplete and code analysis. Recommended for sensitive data +# refer to https://docs.cursor.com/context/ignore-files +.cursorignore +.cursorindexingignore + +# Marimo +marimo/_static/ +marimo/_lsp/ +__marimo__/ diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 0000000..ba1bc74 --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,42 @@ +repos: + - repo: https://github.com/pre-commit/pre-commit-hooks + rev: v4.1.0 + hooks: + - id: double-quote-string-fixer + - id: end-of-file-fixer + exclude: &exclude >- + (?x)^( + tests/.*fixtures| + .*\.json| + .*\.txt + )$ + - id: mixed-line-ending + - id: trailing-whitespace + exclude: *exclude + + - repo: https://github.com/ikamensh/flynt/ + rev: "0.76" + hooks: + - id: flynt + + - repo: https://github.com/pycqa/isort + rev: "5.12.0" + hooks: + - id: isort + exclude: *exclude + + # - repo: local + # hooks: + # - id: pylint + # name: pylint + # entry: pylint + # language: system + # types: [python] + # exclude: > + # (?x)^( + # manage\.py| + # mysite/.*| + # .*/migrations/.*| + # docs/.*| + # tests/.*(? +pip install . +``` ## Description @@ -23,34 +25,69 @@ The simulation models capillary infiltration in porous media using: - **Solver**: Newton-Raphson with adaptive timestepping - **Output**: VTX files for visualization and numpy arrays for data analysis -## Requirements +## External Requirements + +Packages not installable via `pip`: -This project requires FEniCSx and related dependencies with: -- `fenics-dolfinx` - Main FEniCS library -- `petsc4py` - Parallel linear algebra -- `mpi4py` - Message passing interface -- `matplotlib` - Plotting and visualization -- `numpy` - Numerical computing -- `tqdm` - Progress bars +- `fenics-dolfinx==0.9.0` - Main FEniCS library + Python bindings ## Usage -Run the simulation by providing the parameter file: +### Progamatically + +Example. + +```python +from aitw_darcy.data import Params +from aitw_darcy.infiltration import run_simulation + +params = Params() +params.XXX = YYY # Set parameters as needed + +run_simulation(params) +``` + +### CLI + +The installation will make available a `aitw-darcy` command line interface. + +- Run `aitw-darcy --help` to see the available commands. +- Run `aitw-darcy infiltration --help` to see all available options for running the simulation. +- Run `aitw-darcy infiltration INPUT_JSON_FILE` to run the simulation using the parameters in `JSON_FILE`. + +Example with the provided parameter file: ```bash -python infiltration.py data.py +aitw-darcy infiltration example/input.json ``` +#### Tab autocompletion + +Enabling tab autocompletion https://click.palletsprojects.com/en/stable/shell-completion/ + +E.G for `bash` run the command + +```bash +eval "$(_AITW_DARCY_MS_COMPLETE=bash_source aitw-darcy)" +``` + +You can also add it to either `~/.bashrc` or, if you are using a virtual environment, to `bin/activate` of the virtual environment to avoid running the command for every new shell. + ## Output -The simulation generates: +The simulation generates the following output files in the specified output directory: + - `{label}_solution.npy` - Complete solution array for analysis - `{label}_output.bp` - VTX mesh files for ParaView visualization - `{label}.log` - Simulation log file +- `reprod.json` - Copy of the input parameters for reproducibility + +where `{label}` is defined in the input parameters. ## Parameters -Key simulation parameters (configurable in `data.py`): +Key simulation parameters: + - Mesh dimensions and resolution - Material properties (porosity, permeability, surface tension) - Fluid properties (viscosity, contact angle) diff --git a/aitw_darcy/__init__.py b/aitw_darcy/__init__.py new file mode 100644 index 0000000..d8a1d74 --- /dev/null +++ b/aitw_darcy/__init__.py @@ -0,0 +1,2 @@ +"""AITW Darcy Module TODO: improve description""" +__version__ = '0.1.0' diff --git a/aitw_darcy/cli/__init__.py b/aitw_darcy/cli/__init__.py new file mode 100644 index 0000000..10fca5d --- /dev/null +++ b/aitw_darcy/cli/__init__.py @@ -0,0 +1,28 @@ +import click as original_click +import rich_click as click +from rich.traceback import install + +install( + suppress=[original_click, click], +) + +@click.group('aitw-darcy') +def cmd_root(): + """ Command line interface for AITW Darcy simulations """ + pass + +@cmd_root.command() +@click.argument('json_file', required=True, type=click.Path(exists=True)) +@click.option('--output_dir', type=click.Path(), help='Output directory') +def infiltration( + json_file: str, + output_dir: str = '.', + ): + """ Run the infiltration simulation command line interface """ + from aitw_darcy.data import Params + from aitw_darcy.infiltration import run_simulation + + params = Params.from_json(json_file) + if output_dir is not None: + params.output_dir = output_dir + run_simulation(params) diff --git a/aitw_darcy/data.py b/aitw_darcy/data.py new file mode 100644 index 0000000..8493c0a --- /dev/null +++ b/aitw_darcy/data.py @@ -0,0 +1,75 @@ +import json +from dataclasses import dataclass + + +@dataclass +class Params: + # Label for the simulation + label: str = '3d_plate' + output_dir: str = '.' + + # Mesh dimensions (m) + xmin: float = 0 + xmax: float = 5e-2 + ymin: float = 0 + ymax: float = xmax + zmin: float = 0 + zmax: float = 1e-3 + # Number of nodes + nx: int = 40 + ny: int = 160 + nz: int = 40 + + # Timestepping + # Initial time step length (s) + delta_t0: float = 1e-1 + delta_tmin: float = 1e-2 + delta_tmax: float = 10 + # Simulation time (s) + t_end: float = 120 + + # Porosity + porosity: float = 0.6 + + # Surface tension (N/m) + surface_tension: float = 30e-3 + + # Contact angle + contact_angle: float = 10 + + # Mean pore radius (m) + mean_pore_radius: float = 20e-6 + + # Relative standard deviation of pore radius + rstd: float = 1/5 + + # Pore radius cutoff + r_cutoff: float = 2e-6 + + # Permeability + k_long: float = 1e-12 + k_trans: float = 1e-16 + + # Fluid viscosity (Pa s) + viscosity: float = 1e-3 + + # External pressure + u_ext: float = -1e-6 + + def to_json(self, json_file: str): + """Save the parameters to a JSON file""" + data = {k: v for k, v in self.__dict__.items() if not k.startswith('_')} + with open(json_file, 'w') as f: + json.dump(data, f, indent=4) + + @classmethod + def from_json(cls, json_file: str) -> 'Params': + """Create an instance from a JSON file""" + with open(json_file, 'r') as f: + data = json.load(f) + return cls.from_dict(data) + + @classmethod + def from_dict(cls, data: dict) -> 'Params': + """Create an instance from a JSON file""" + return cls(**data) diff --git a/aitw_darcy/infiltration.py b/aitw_darcy/infiltration.py new file mode 100644 index 0000000..5f5fcb3 --- /dev/null +++ b/aitw_darcy/infiltration.py @@ -0,0 +1,249 @@ +"""Module for running infiltration simulation using Darcy's law.""" +import os + +import numpy as np +import ufl +from basix.ufl import element +from dolfinx import default_scalar_type, fem, io, log, mesh +from dolfinx.fem.petsc import NonlinearProblem +from dolfinx.nls.petsc import NewtonSolver +# from matplotlib import pyplot +from mpi4py import MPI +# from petsc4py import PETSc +from tqdm import tqdm +from ufl import cos, dot, dx, erf, grad, pi, sqrt + +from .data import Params + + +def run_simulation(parameters: Params): + """ Run the infiltration simulation with given parameters""" + label = parameters.label + xmin = parameters.xmin + xmax = parameters.xmax + ymin = parameters.ymin + ymax = parameters.ymax + zmin = parameters.zmin + zmax = parameters.zmax + nx = parameters.nx + ny = parameters.ny + nz = parameters.nz + delta_t0 = parameters.delta_t0 + delta_tmin = parameters.delta_tmin + delta_tmax = parameters.delta_tmax + t_end = parameters.t_end + + porosity_phi = parameters.porosity + surf_tens_gamma = parameters.surface_tension + c_ang = parameters.contact_angle + + mean_pore_radius = parameters.mean_pore_radius + rstd = parameters.rstd + + r_cutoff = parameters.r_cutoff + + k_long = parameters.k_long + k_trans = parameters.k_trans + + viscosity = parameters.viscosity + u_ext = parameters.u_ext + + output_dir = parameters.output_dir + + # Output + solution_file = f'{label}_solution.npy' # Full solution + output_file = f'{label}_output.bp' # VTX mesh for visualization + + # Ensure output directory exists + os.makedirs(output_dir, exist_ok=True) + solution_file = os.path.join(output_dir, solution_file) + output_file = os.path.join(output_dir, output_file) + log_file = os.path.join(output_dir, f'{label}.log') + + # Save input parameters to output directory for reproducibility + reprod_file = os.path.join(output_dir, f'reprod.json') + parameters.to_json(reprod_file) + + log.set_output_file(log_file) + log.set_log_level(log.LogLevel.INFO) + + # Helper functions + + def get_variable_value(v): + """ Get the nodal values for an expression v """ + q = fem.Function(V) + q.interpolate(fem.Expression(v, V.element.interpolation_points())) + return q.x.array + + def ramp_up(val=1, t1=0.1): + """ Linear ramp from zero beginning for t = 0, useful to assist convergence in the beginning """ + return val*ufl.min_value(1, t/t1) + + # Generate 3D mesh + domain = mesh.create_box( + MPI.COMM_WORLD, + [[xmin, ymin, zmin], [xmax, ymax, zmax]], [nx, ny, nz], + mesh.CellType.hexahedron + ) + + x = ufl.SpatialCoordinate(domain) + + # Current time and time step size + t = fem.Constant(domain, default_scalar_type(0)) + delta_t = fem.Constant(domain, default_scalar_type(delta_t0)) + + # Create function space + P1 = element('Lagrange', domain.basix_cell(), 1) + V = fem.functionspace(domain, P1) + + # Unknown solution (t=n+1) + u = fem.Function(V, name='u') + + # Previous time step solution + u_n = fem.Function(V, name='u_n') + + # Saturation curve + def S_curve(r): + sigma = rstd*mean_pore_radius + return 1/2*(erf((r-mean_pore_radius)/(sqrt(2)*sigma)) - erf((r_cutoff-mean_pore_radius)/(sqrt(2)*sigma))) + + # Capillary radius + theta = c_ang*pi/180 + r_c = -2*surf_tens_gamma*cos(theta)/u + + # Saturation + #S = ufl.conditional(u < 0, S_curve(r_c), 1.0) + S = S_curve(r_c) + + # Previous saturation + S_n = ufl.replace(S, {u: u_n}) + + # Flux + k = ufl.as_tensor([[k_long, 0, 0], [0, k_trans, 0], [0, 0, k_trans]]) + k_s = 1 + ramp_up(-1+S) + J = -k_s*k/viscosity*grad(u) + + ############### + # Weak form + ############### + + delta_u = ufl.TestFunction(V) + + F = delta_u * porosity_phi*(S - S_n) * dx - delta_t*dot(grad(delta_u), J) * dx + + # Direchlet boundary conditions + + fdim = domain.topology.dim - 1 + domain.topology.create_connectivity(fdim, fdim+1) + boundary_facets = mesh.exterior_facet_indices(domain.topology) + boundary_dofs = fem.locate_dofs_topological(V, fdim, boundary_facets) + bc = fem.dirichletbc(default_scalar_type(u_ext), boundary_dofs, V) + bcs = [bc] + + # Initial pressure + u0 = -2*surf_tens_gamma*cos(theta)/r_cutoff + + # Initial moisture content + u.interpolate(lambda x: np.full(x[0].shape, u0)) + u.x.array[boundary_dofs] = u_ext + + # Initialize the Newton solver + problem = NonlinearProblem(F, u, bcs) + # solver = NewtonSolver(MPI.COMM_WORLD, problem) + solver = NewtonSolver(MPI.COMM_WORLD, problem) + # Solution is rather sensitive to tolerance, + # numerical accuracy seems to be close to 1e-11 + solver.rtol = 0 + solver.atol = 1e-11 + solver.report = True + solver.error_on_nonconvergence = False + solver.max_it = 10 + solver.relaxation_parameter = 0.8 + + # Use mumps LU solver + ksp = solver.krylov_solver + ksp.setType('preonly') + ksp.getPC().setType('lu') + ksp.getPC().setFactorSolverType('mumps') + + # Solutions at each time step + sol = [u.x.array.copy()] + + # Time values + times = [t.value.item()] + + if domain.comm.rank == 0: + pbar = tqdm(total=t_end, unit_scale=True) + + S_func = fem.Function(V, name='S') + + with io.VTXWriter(domain.comm, output_file, [u, S_func]) as vtxfile: + while True: + u_n.x.array[:] = u.x.array + + # Adaptive timestepping + while True: + n, converged = solver.solve(u) + if converged: + solver.relaxation_parameter = 0.8 + if n <= 3: + delta_t.value = min(delta_t.value * 1.2, delta_tmax) + log.log(log.LogLevel.INFO, f'Timestep = {delta_t.value}') + elif n >= 7: + delta_t.value = delta_t.value * 0.5 + log.log(log.LogLevel.INFO, f'Timestep = {delta_t.value}') + break + else: + solver.relaxation_parameter = 0.4 + delta_t.value = delta_t.value * 0.5 + if delta_t.value < delta_tmin: + raise RuntimeError('Minimum timestep size reached!') + u.x.array[:] = u_n.x.array + log.log(log.LogLevel.INFO, f'Timestep = {delta_t.value}') + + t.value += delta_t.value + + # Record values + sol.append(u.x.array.copy()) + times.append(t.value.item()) + + if domain.comm.rank == 0: + pbar.update(delta_t.value) + + # Save output + S_func.interpolate(fem.Expression(S, V.element.interpolation_points())) + vtxfile.write(t.value.item()) + + # Termination criterion + if t.value > t_end: + break + + S_threshold = 0.99 + S_vals = get_variable_value(S) + finished = np.all(S_vals > S_threshold) + statuses = domain.comm.gather(finished, root=0) + can_stop = None + if domain.comm.rank == 0: + can_stop = all(statuses) + if domain.comm.bcast(can_stop, root=0): + break + + # Numpy array output for single parallel runs + if domain.comm.size == 1: + sol = np.array(sol) + np.save(solution_file, sol) + + +# The residual and Jacobian for diagnostics +#residual = fem.form(F) +#RD = fem.petsc.create_vector(residual) +#fem.petsc.assemble_vector(RD, residual) +#jacobian = fem.form(ufl.derivative(F, u)) + +#fem.petsc.apply_lifting(RD, [jacobian], [bcs]) + +#A = fem.petsc.create_matrix(jacobian) +#fem.petsc.assemble_matrix(A, jacobian) +#A.assemble() +#print(domain.comm.rank, A.getSize()) +#m = A.convert('dense').getDenseArray() diff --git a/data.py b/data.py deleted file mode 100644 index e4a8430..0000000 --- a/data.py +++ /dev/null @@ -1,57 +0,0 @@ -### -# -# Input parameters -# -### - -# Label for the simulation -label = '3d_plate' - -# Mesh dimensions (m) -xmin = 0 -xmax = 5e-2 -ymin = 0 -ymax = xmax -zmin = 0 -zmax = 1e-3 -nx = 40 # Number of nodes -ny = 160 -nz = 40 - -# Timestepping -Δt0 = 1e-1 # Initial time step length (s) -Δt_min = 1e-2 -Δt_max = 10 -t_end = 120 # Simulation time (s) - -# Porosity -ϕ = 0.6 - -# Surface tension (N/m) -γ = 30e-3 - -# Contact angle -θdeg = 10 - -# Mean pore radius (m) -r_μ = 20e-6 - -# Relative standard deviation of pore radius -rstd = 1/5 - -# Pore radius cutoff -r_cutoff = 2e-6 - - -# Permeability -k_long = 1e-12 -k_trans = 1e-16 - -# Fluid viscosity (Pa s) -μ = 1e-3 - -# External pressure -u_ext = -1e-6 - -### -### diff --git a/example/input.json b/example/input.json new file mode 100644 index 0000000..a288fe8 --- /dev/null +++ b/example/input.json @@ -0,0 +1,26 @@ +{ + "label": "3d_plate", + "xmin": 0, + "xmax": 0.05, + "ymin": 0, + "ymax": 0.05, + "zmin": 0, + "zmax": 0.001, + "nx": 40, + "ny": 160, + "nz": 40, + "delta_t0": 0.1, + "delta_tmin": 0.01, + "delta_tmax": 10, + "t_end": 120, + "porosity": 0.6, + "surface_tension": 0.03, + "contact_angle": 10, + "mean_pore_radius": 2e-05, + "rstd": 0.2, + "r_cutoff": 2e-06, + "k_long": 1e-12, + "k_trans": 1e-16, + "viscosity": 0.001, + "u_ext": -1e-06 +} \ No newline at end of file diff --git a/infiltration.py b/infiltration.py deleted file mode 100644 index 57f9f25..0000000 --- a/infiltration.py +++ /dev/null @@ -1,214 +0,0 @@ -from matplotlib import pyplot -import numpy as np -from tqdm import trange, tqdm - -from petsc4py import PETSc -from mpi4py import MPI - -import ufl -from ufl import dot, grad, dx, exp, inner, sqrt, erf, cos, pi -from dolfinx import fem, mesh, io, default_scalar_type, log -from dolfinx.fem.petsc import NonlinearProblem -from dolfinx.nls.petsc import NewtonSolver -from basix.ufl import element - -import sys -from runpy import run_path - -# Inject the variables from the parameter file to the global namespace -try: - param_file = sys.argv[1] - params = run_path(param_file) -except: - error("Please provide the parameter file (.py) in the working directory!") - -for k, v in params.items(): - if not k.startswith('_'): - globals()[k] = v - - -# Output -solution_file = f'{label}_solution.npy' # Full solution -output_file = f'{label}_output.bp' # VTX mesh for visualization - -log.set_output_file(f'{label}.log') -log.set_log_level(log.LogLevel.INFO) - -# Helper functions - -def get_variable_value(v): - """ Get the nodal values for an expression v """ - q = fem.Function(V) - q.interpolate(fem.Expression(v, V.element.interpolation_points())) - return q.x.array - -def ramp_up(val=1, t1=0.1): - """ Linear ramp from zero beginning for t = 0, useful to assist convergence in the beginning """ - return val*ufl.min_value(1, t/t1) - -# Generate 3D mesh -domain = mesh.create_box(MPI.COMM_WORLD, [[xmin, ymin, zmin], [xmax, ymax, zmax]], [nx, ny, nz], - mesh.CellType.hexahedron) - -x = ufl.SpatialCoordinate(domain) - -# Current time and time step size -t = fem.Constant(domain, default_scalar_type(0)) -Δt = fem.Constant(domain, default_scalar_type(Δt0)) - -# Create function space -P1 = element("Lagrange", domain.basix_cell(), 1) -V = fem.functionspace(domain, P1) - -# Unknown solution (t=n+1) -u = fem.Function(V, name='u') - -# Previous time step solution -u_n = fem.Function(V, name='u_n') - -# Saturation curve -def S_curve(r): - σ = rstd*r_μ - return 1/2*(erf((r-r_μ)/(sqrt(2)*σ)) - erf((r_cutoff-r_μ)/(sqrt(2)*σ))) - -# Capillary radius -θ = θdeg*pi/180 -r_c = -2*γ*cos(θ)/u - -# Saturation -#S = ufl.conditional(u < 0, S_curve(r_c), 1.0) -S = S_curve(r_c) - -# Previous saturation -S_n = ufl.replace(S, {u: u_n}) - -# Flux -k = ufl.as_tensor([[k_long, 0, 0], [0, k_trans, 0], [0, 0, k_trans]]) -k_s = 1 + ramp_up(-1+S) -J = -k_s*k/μ*grad(u) - -############### -# Weak form -############### - -δu = ufl.TestFunction(V) - -F = δu * ϕ*(S - S_n) * dx - Δt*dot(grad(δu), J) * dx - -# Direchlet boundary conditions - -fdim = domain.topology.dim - 1 -domain.topology.create_connectivity(fdim, fdim+1) -boundary_facets = mesh.exterior_facet_indices(domain.topology) -boundary_dofs = fem.locate_dofs_topological(V, fdim, boundary_facets) -bc = fem.dirichletbc(default_scalar_type(u_ext), boundary_dofs, V) -bcs = [bc] - -# Initial pressure -u0 = -2*γ*cos(θ)/r_cutoff - -if __name__ == '__main__': - # Initial moisture content - u.interpolate(lambda x: np.full(x[0].shape, u0)) - u.x.array[boundary_dofs] = u_ext - - # Initialize the Newton solver - problem = NonlinearProblem(F, u, bcs) - solver = NewtonSolver(MPI.COMM_WORLD, problem) - - # Solution is rather sensitive to tolerance, - # numerical accuracy seems to be close to 1e-11 - solver.rtol = 0 - solver.atol = 1e-11 - solver.report = True - solver.error_on_nonconvergence = False - solver.max_it = 10 - solver.relaxation_parameter = 0.8 - - # Use mumps LU solver - ksp = solver.krylov_solver - ksp.setType("preonly") - ksp.getPC().setType("lu") - ksp.getPC().setFactorSolverType("mumps") - - # Solutions at each time step - sol = [u.x.array.copy()] - - # Time values - times = [t.value.item()] - - if domain.comm.rank == 0: - pbar = tqdm(total=t_end, unit_scale=True) - - S_func = fem.Function(V, name='S') - - with io.VTXWriter(domain.comm, output_file, [u, S_func]) as vtxfile: - while True: - u_n.x.array[:] = u.x.array - - # Adaptive timestepping - while True: - n, converged = solver.solve(u) - if converged: - solver.relaxation_parameter = 0.8 - if n <= 3: - Δt.value = min(Δt.value * 1.2, Δt_max) - log.log(log.LogLevel.INFO, f'Timestep = {Δt.value}') - elif n >= 7: - Δt.value = Δt.value * 0.5 - log.log(log.LogLevel.INFO, f'Timestep = {Δt.value}') - break - else: - solver.relaxation_parameter = 0.4 - Δt.value = Δt.value * 0.5 - if Δt.value < Δt_min: - raise Error("Minimum timestep size reached!") - u.x.array[:] = u_n.x.array - log.log(log.LogLevel.INFO, f'Timestep = {Δt.value}') - - t.value += Δt.value - - # Record values - sol.append(u.x.array.copy()) - times.append(t.value.item()) - - if domain.comm.rank == 0: - pbar.update(Δt.value) - - # Save output - S_func.interpolate(fem.Expression(S, V.element.interpolation_points())) - vtxfile.write(t.value.item()) - - # Termination criterion - if t.value > t_end: - break - - S_threshold = 0.99 - S_vals = get_variable_value(S) - finished = np.all(S_vals > S_threshold) - statuses = domain.comm.gather(finished, root=0) - can_stop = None - if domain.comm.rank == 0: - can_stop = all(statuses) - if domain.comm.bcast(can_stop, root=0): - break - - # Numpy array output for single parallel runs - if domain.comm.size == 1: - sol = np.array(sol) - np.save(solution_file, sol) - - -# The residual and Jacobian for diagnostics -#residual = fem.form(F) -#RD = fem.petsc.create_vector(residual) -#fem.petsc.assemble_vector(RD, residual) -#jacobian = fem.form(ufl.derivative(F, u)) - -#fem.petsc.apply_lifting(RD, [jacobian], [bcs]) - -#A = fem.petsc.create_matrix(jacobian) -#fem.petsc.assemble_matrix(A, jacobian) -#A.assemble() -#print(domain.comm.rank, A.getSize()) -#m = A.convert('dense').getDenseArray() diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..ae6e564 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,56 @@ +[build-system] +# build the package with [flit](https://flit.readthedocs.io) +requires = ["flit_core >=3.4,<4"] +build-backend = "flit_core.buildapi" + +[project] +# See https://www.python.org/dev/peps/pep-0621/ +name = "AITW_darcy_infiltration" +dynamic = ["version"] +description = "Code for simulating infiltration processes in porous materials using the FEniCSx finite element framework." +authors = [ + { name = "Royson Donate Dsouza" } +] +readme = "README.md" +license = { file = "LICENSE" } +classifiers = [ + "Development Status :: 2 - Pre-Alpha", + "Operating System :: Microsoft :: Windows", + "Operating System :: POSIX :: Linux", + "Operating System :: MacOS", + "Programming Language :: Python", +] +keywords = ["finite-elements", "darcy", "infiltration", "multiscale"] +requires-python = ">=3.10" +dependencies = [ + "rich_click", + "tqdm", + "numpy", + + "mpi4py", + "petsc4py", + "fenics-ufl", + "fenics-basix==0.9.0", +] + +[project.urls] +Source = "https://github.com/VTT-ProperTune/AITW-Darcy-Infiltration-Model" + +[project.optional-dependencies] +pre-commit = [ + "pre-commit", +] +release = [ + "flit" +] + +[project.scripts] +aitw-darcy = "aitw_darcy.cli:cmd_root" + +[tool.flit.module] +name = "aitw_darcy" + +[tool.flit.sdist] +exclude = [ + ".gitignore", ".github", ".pre-commit-config.yaml", +]