diff --git a/.devcontainer/Dockerfile b/.devcontainer/Dockerfile index 790ed16..1ff81d1 100644 --- a/.devcontainer/Dockerfile +++ b/.devcontainer/Dockerfile @@ -1,4 +1,4 @@ -FROM dolfinx/dolfinx:v0.7.0 +FROM dolfinx/dolfinx:v0.6.0 # FROM mcr.microsoft.com/devcontainers/base:jammy ARG DEBIAN_FRONTEND=noninteractive # RUN useradd -ms /bin/bash vscode @@ -25,23 +25,22 @@ RUN DEBIAN_FRONTEND=noninteractive \ libffi-dev \ liblzma-dev -# Python and poetry installation -USER $USER -ARG HOME="/home/$USER" +# # Python and poetry installation +# # USER $USER +# # ARG HOME="/home/$USER" -# Don't change. Dolfinx in the base image is only compiled with 3.10 -ARG PYTHON_VERSION=3.10 - -# Use poetry pyenv -ENV PYENV_ROOT="${HOME}/.pyenv" -ENV PATH="${PYENV_ROOT}/shims:${PYENV_ROOT}/bin:${HOME}/.local/bin:$PATH" +# # Don't change. Dolfinx in the base image is only compiled with 3.10 +# #ARG PYTHON_VERSION=3.10.6 +# # Use poetry pyenv +# #ENV PYENV_ROOT="${HOME}/.pyenv" +# #ENV PATH="${PYENV_ROOT}/shims:${PYENV_ROOT}/bin:${HOME}/.local/bin:/root/.local/bin:$PATH" +ENV PATH="/root/.local/bin:$PATH" RUN echo "done 0" \ - && curl https://pyenv.run | bash \ - && echo "done 1" \ - && pyenv install ${PYTHON_VERSION} \ - && echo "done 2" \ - && pyenv global ${PYTHON_VERSION} \ - && echo "done 3" \ - && curl -sSL https://install.python-poetry.org | python3 - \ - && poetry config virtualenvs.in-project true +# # && curl https://pyenv.run | bash \ +# # && echo "done 1" \ +# # && pyenv install ${PYTHON_VERSION} \ +# # && echo "done 2" \ +# # && pyenv global ${PYTHON_VERSION} \ +# # && echo "done 3" \ + && curl -sSL https://install.python-poetry.org | python3 - diff --git a/.devcontainer/postCreate.sh b/.devcontainer/postCreate.sh index af415a5..3e4b744 100644 --- a/.devcontainer/postCreate.sh +++ b/.devcontainer/postCreate.sh @@ -1,2 +1,4 @@ +poetry config virtualenvs.create false +poetry env use system poetry install --with dev --no-interaction -poetry run pip install petsc4py==3.20.0 mpi4py==3.1.5 https://github.com/FEniCS/ufl/archive/refs/tags/2023.2.0.zip https://github.com/FEniCS/ffcx/archive/v0.7.0.zip https://github.com/FEniCS/basix/archive/v0.7.0.zip \ No newline at end of file +#poetry run pip install petsc4py==3.20.0 mpi4py==3.1.5 https://github.com/FEniCS/ufl/archive/refs/tags/2023.2.0.zip https://github.com/FEniCS/ffcx/archive/v0.7.0.zip https://github.com/FEniCS/basix/archive/v0.7.0.zip \ No newline at end of file diff --git a/.github/workflows/python-test-3d.yml b/.github/workflows/python-test-3d.yml index 53d758d..5e1d915 100644 --- a/.github/workflows/python-test-3d.yml +++ b/.github/workflows/python-test-3d.yml @@ -4,6 +4,7 @@ name: Tests3D on: + workflow_dispatch: workflow_call: inputs: event_type: @@ -23,12 +24,11 @@ jobs: - name: Run uses: tj-actions/docker-run@v2 with: - image: dolfinx/dolfinx:v0.7.0 + image: dolfinx/dolfinx:v0.6.0 name: dolfinx options: -v ${{ github.workspace }}:/home/warmth args: | - bash -c "cd /home/warmth" && pip install . pytest==7.4.2 pytest-cov==4.1.0 && pytest --cov-report=term-missing --cov=warmth/3d tests/3d | tee pytest-coverage.txt - + bash -c "cd /home/warmth" && pip install . pytest==7.4.2 pytest-cov==4.1.0 && pytest --cov-report=term-missing --cov=warmth tests | tee pytest-coverage.txt - name: Comment coverage if: ${{ github.event_name == 'pull_request' && github.event.action == 'opened' }} diff --git a/docs/notebooks/data/0.gri b/docs/notebooks/data/0.gri index 120ca3b..37d352e 100644 Binary files a/docs/notebooks/data/0.gri and b/docs/notebooks/data/0.gri differ diff --git a/docs/notebooks/data/100.gri b/docs/notebooks/data/100.gri index b899d47..1b24b51 100644 Binary files a/docs/notebooks/data/100.gri and b/docs/notebooks/data/100.gri differ diff --git a/docs/notebooks/data/163.gri b/docs/notebooks/data/163.gri index 8bbc9a4..adab1d6 100644 Binary files a/docs/notebooks/data/163.gri and b/docs/notebooks/data/163.gri differ diff --git a/docs/notebooks/data/168.gri b/docs/notebooks/data/168.gri index 680f5c4..64f91b8 100644 Binary files a/docs/notebooks/data/168.gri and b/docs/notebooks/data/168.gri differ diff --git a/docs/notebooks/data/170.gri b/docs/notebooks/data/170.gri index 9e6d4aa..90529c4 100644 Binary files a/docs/notebooks/data/170.gri and b/docs/notebooks/data/170.gri differ diff --git a/docs/notebooks/data/182.gri b/docs/notebooks/data/182.gri index cbd108f..496b294 100644 Binary files a/docs/notebooks/data/182.gri and b/docs/notebooks/data/182.gri differ diff --git a/docs/notebooks/data/66.gri b/docs/notebooks/data/66.gri index a075e84..4f0a8a3 100644 Binary files a/docs/notebooks/data/66.gri and b/docs/notebooks/data/66.gri differ diff --git a/tests/__test_dol.py b/tests/__test_dol.py deleted file mode 100644 index 015a761..0000000 --- a/tests/__test_dol.py +++ /dev/null @@ -1,322 +0,0 @@ -# Copyright (C) 2019 Chris Richardson -# -# This file is part of DOLFINx (https://www.fenicsproject.org) -# -# SPDX-License-Identifier: LGPL-3.0-or-later -"""Unit tests for assembly over domains""" - -from mpi4py import MPI - -import numpy as np -import pytest - -import ufl -from dolfinx import cpp as _cpp -from dolfinx import default_scalar_type, fem, la -from dolfinx.fem import (Constant, Function, assemble_scalar, dirichletbc, - form, functionspace) -from dolfinx.mesh import (GhostMode, Mesh, create_unit_square, locate_entities, - locate_entities_boundary, meshtags, - meshtags_from_entities) - - -@pytest.fixture -def mesh(): - return create_unit_square(MPI.COMM_WORLD, 10, 10) - - -def create_cell_meshtags_from_entities(mesh: Mesh, dim: int, cells: np.ndarray, values: np.ndarray): - mesh.topology.create_connectivity(mesh.topology.dim, 0) - cell_to_vertices = mesh.topology.connectivity(mesh.topology.dim, 0) - entities = _cpp.graph.AdjacencyList_int32([cell_to_vertices.links(cell) for cell in cells]) - return meshtags_from_entities(mesh, dim, entities, values) - - -parametrize_ghost_mode = pytest.mark.parametrize("mode", [ - pytest.param(GhostMode.none, marks=pytest.mark.skipif(condition=MPI.COMM_WORLD.size > 1, - reason="Unghosted interior facets fail in parallel")), - pytest.param(GhostMode.shared_facet, marks=pytest.mark.skipif(condition=MPI.COMM_WORLD.size == 1, - reason="Shared ghost modes fail in serial"))]) - - -@pytest.mark.parametrize("mode", [GhostMode.none, GhostMode.shared_facet]) -@pytest.mark.parametrize("meshtags_factory", [meshtags, create_cell_meshtags_from_entities]) -def test_assembly_dx_domains(mode, meshtags_factory): - mesh = create_unit_square(MPI.COMM_WORLD, 10, 10, ghost_mode=mode) - V = functionspace(mesh, ("Lagrange", 1)) - u, v = ufl.TrialFunction(V), ufl.TestFunction(V) - - # Prepare a marking structures - # indices cover all cells - # values are [1, 2, 3, 3, ...] - cell_map = mesh.topology.index_map(mesh.topology.dim) - num_cells = cell_map.size_local + cell_map.num_ghosts - indices = np.arange(0, num_cells) - values = np.full(indices.shape, 3, dtype=np.intc) - values[0] = 1 - values[1] = 2 - marker = meshtags_factory(mesh, mesh.topology.dim, indices, values) - dx = ufl.Measure('dx', subdomain_data=marker, domain=mesh) - w = Function(V) - w.x.array[:] = 0.5 - - # Assemble matrix - a = form(w * ufl.inner(u, v) * (dx(1) + dx(2) + dx(3))) - A = fem.assemble_matrix(a) - A.scatter_reverse() - a2 = form(w * ufl.inner(u, v) * dx) - A2 = fem.assemble_matrix(a2) - A2.scatter_reverse() - assert np.allclose(A.data, A2.data) - - bc = dirichletbc(Function(V), range(30)) - - # Assemble vector - L = form(ufl.inner(w, v) * (dx(1) + dx(2) + dx(3))) - b = fem.assemble_vector(L) - - fem.apply_lifting(b.array, [a], [[bc]]) - b.scatter_reverse(la.InsertMode.add) - fem.set_bc(b.array, [bc]) - - L2 = form(ufl.inner(w, v) * dx) - b2 = fem.assemble_vector(L2) - fem.apply_lifting(b2.array, [a], [[bc]]) - b2.scatter_reverse(la.InsertMode.add) - fem.set_bc(b2.array, [bc]) - assert np.allclose(b.array, b2.array) - - # Assemble scalar - L = form(w * (dx(1) + dx(2) + dx(3))) - s = assemble_scalar(L) - s = mesh.comm.allreduce(s, op=MPI.SUM) - assert s == pytest.approx(0.5, rel=1.0e-6) - L2 = form(w * dx) - s2 = assemble_scalar(L2) - s2 = mesh.comm.allreduce(s2, op=MPI.SUM) - assert s == pytest.approx(s2, rel=1.0e-6) - - -@pytest.mark.parametrize("mode", [GhostMode.none, GhostMode.shared_facet]) -def test_assembly_ds_domains(mode): - mesh = create_unit_square(MPI.COMM_WORLD, 10, 10, ghost_mode=mode) - V = functionspace(mesh, ("Lagrange", 1)) - u, v = ufl.TrialFunction(V), ufl.TestFunction(V) - - def bottom(x): - return np.isclose(x[1], 0.0) - - def top(x): - return np.isclose(x[1], 1.0) - - def left(x): - return np.isclose(x[0], 0.0) - - def right(x): - return np.isclose(x[0], 1.0) - - bottom_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, bottom) - bottom_vals = np.full(bottom_facets.shape, 1, np.intc) - - top_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, top) - top_vals = np.full(top_facets.shape, 2, np.intc) - - left_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, left) - left_vals = np.full(left_facets.shape, 3, np.intc) - - right_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, right) - right_vals = np.full(right_facets.shape, 6, np.intc) - - indices = np.hstack((bottom_facets, top_facets, left_facets, right_facets)) - values = np.hstack((bottom_vals, top_vals, left_vals, right_vals)) - - indices, pos = np.unique(indices, return_index=True) - marker = meshtags(mesh, mesh.topology.dim - 1, indices, values[pos]) - - ds = ufl.Measure('ds', subdomain_data=marker, domain=mesh) - - w = Function(V) - w.x.array[:] = 0.5 - - bc = dirichletbc(Function(V), range(30)) - - # Assemble matrix - a = form(w * ufl.inner(u, v) * (ds(1) + ds(2) + ds(3) + ds(6))) - A = fem.assemble_matrix(a) - A.scatter_reverse() - a2 = form(w * ufl.inner(u, v) * ds) - A2 = fem.assemble_matrix(a2) - A2.scatter_reverse() - assert np.allclose(A.data, A2.data) - - # Assemble vector - L = form(ufl.inner(w, v) * (ds(1) + ds(2) + ds(3) + ds(6))) - b = fem.assemble_vector(L) - - fem.apply_lifting(b.array, [a], [[bc]]) - b.scatter_reverse(la.InsertMode.add) - fem.set_bc(b.array, [bc]) - - L2 = form(ufl.inner(w, v) * ds) - b2 = fem.assemble_vector(L2) - fem.apply_lifting(b2.array, [a2], [[bc]]) - b2.scatter_reverse(la.InsertMode.add) - fem.set_bc(b2.array, [bc]) - assert np.allclose(b.array, b2.array) - - # Assemble scalar - L = form(w * (ds(1) + ds(2) + ds(3) + ds(6))) - s = assemble_scalar(L) - s = mesh.comm.allreduce(s, op=MPI.SUM) - L2 = form(w * ds) - s2 = assemble_scalar(L2) - s2 = mesh.comm.allreduce(s2, op=MPI.SUM) - assert s == pytest.approx(s2, 1.0e-6) - assert 2.0 == pytest.approx(s, 1.0e-6) # /NOSONAR - - -@parametrize_ghost_mode -def test_assembly_dS_domains(mode): - N = 10 - mesh = create_unit_square(MPI.COMM_WORLD, N, N, ghost_mode=mode) - one = Constant(mesh, default_scalar_type(1)) - val = assemble_scalar(form(one * ufl.dS)) - val = mesh.comm.allreduce(val, op=MPI.SUM) - assert val == pytest.approx(2 * (N - 1) + N * np.sqrt(2), 1.0e-5) - - -@parametrize_ghost_mode -def test_additivity(mode): - mesh = create_unit_square(MPI.COMM_WORLD, 12, 12, ghost_mode=mode) - V = functionspace(mesh, ("Lagrange", 1)) - - f1 = Function(V) - f2 = Function(V) - f3 = Function(V) - f1.x.array[:] = 1.0 - f2.x.array[:] = 2.0 - f3.x.array[:] = 3.0 - j1 = ufl.inner(f1, f1) * ufl.dx(mesh) - j2 = ufl.inner(f2, f2) * ufl.ds(mesh) - j3 = ufl.inner(ufl.avg(f3), ufl.avg(f3)) * ufl.dS(mesh) - - # Assemble each scalar form separately - J1 = mesh.comm.allreduce(assemble_scalar(form(j1)), op=MPI.SUM) - J2 = mesh.comm.allreduce(assemble_scalar(form(j2)), op=MPI.SUM) - J3 = mesh.comm.allreduce(assemble_scalar(form(j3)), op=MPI.SUM) - - # Sum forms and assemble the result - J12 = mesh.comm.allreduce(assemble_scalar(form(j1 + j2)), op=MPI.SUM) - J13 = mesh.comm.allreduce(assemble_scalar(form(j1 + j3)), op=MPI.SUM) - J23 = mesh.comm.allreduce(assemble_scalar(form(j2 + j3)), op=MPI.SUM) - J123 = mesh.comm.allreduce(assemble_scalar(form(j1 + j2 + j3)), op=MPI.SUM) - - # Compare assembled values - assert (J1 + J2) == pytest.approx(J12) - assert (J1 + J3) == pytest.approx(J13) - assert (J2 + J3) == pytest.approx(J23) - assert (J1 + J2 + J3) == pytest.approx(J123) - - -def test_manual_integration_domains(): - """Test that specifying integration domains manually i.e. - by passing a list of cell indices or (cell, local facet) pairs to - form gives the same result as the usual approach of tagging""" - n = 4 - msh = create_unit_square(MPI.COMM_WORLD, n, n) - - V = functionspace(msh, ("Lagrange", 1)) - u = ufl.TrialFunction(V) - v = ufl.TestFunction(V) - - # Create meshtags to mark some cells - tdim = msh.topology.dim - cell_map = msh.topology.index_map(tdim) - num_cells = cell_map.size_local + cell_map.num_ghosts - cell_indices = np.arange(0, num_cells) - cell_values = np.zeros_like(cell_indices, dtype=np.intc) - marked_cells = locate_entities(msh, tdim, lambda x: x[0] < 0.75) - cell_values[marked_cells] = 7 - mt_cells = meshtags(msh, tdim, cell_indices, cell_values) - - # Create meshtags to mark some exterior facets - msh.topology.create_entities(tdim - 1) - facet_map = msh.topology.index_map(tdim - 1) - num_facets = facet_map.size_local + facet_map.num_ghosts - facet_indices = np.arange(0, num_facets) - facet_values = np.zeros_like(facet_indices, dtype=np.intc) - marked_ext_facets = locate_entities_boundary(msh, tdim - 1, lambda x: np.isclose(x[0], 0.0)) - marked_int_facets = locate_entities(msh, tdim - 1, lambda x: x[0] < 0.75) - # marked_int_facets will also contain facets on the boundary, so set - # these values first, followed by the values for marked_ext_facets - facet_values[marked_int_facets] = 3 - facet_values[marked_ext_facets] = 6 - mt_facets = meshtags(msh, tdim - 1, facet_indices, facet_values) - - # Create measures - dx_mt = ufl.Measure("dx", subdomain_data=mt_cells, domain=msh) - ds_mt = ufl.Measure("ds", subdomain_data=mt_facets, domain=msh) - dS_mt = ufl.Measure("dS", subdomain_data=mt_facets, domain=msh) - - g = Function(V) - g.interpolate(lambda x: x[1]**2) - - def create_forms(dx, ds, dS): - a = form(ufl.inner(g * u, v) * (dx(0) + dx(7) + ds(6)) + ufl.inner(g * u("+"), v("+") + v("-")) * dS(3)) - L = form(ufl.inner(g, v) * (dx(0) + dx(7) + ds(6)) + ufl.inner(g, v("+") + v("-")) * dS(3)) - return (a, L) - - # Create forms and assemble - a, L = create_forms(dx_mt, ds_mt, dS_mt) - A_mt = fem.assemble_matrix(a) - A_mt.scatter_reverse() - b_mt = fem.assemble_vector(L) - - # Manually specify cells to integrate over (removing ghosts - # to give same result as above) - cell_domains = [(domain_id, cell_indices[(cell_values == domain_id) - & (cell_indices < cell_map.size_local)]) - for domain_id in [0, 7]] - - # Manually specify exterior facets to integrate over as - # (cell, local facet) pairs - ext_facet_domain = [] - msh.topology.create_connectivity(tdim, tdim - 1) - msh.topology.create_connectivity(tdim - 1, tdim) - c_to_f = msh.topology.connectivity(tdim, tdim - 1) - f_to_c = msh.topology.connectivity(tdim - 1, tdim) - for f in marked_ext_facets: - if f < facet_map.size_local: - c = f_to_c.links(f)[0] - local_f = np.where(c_to_f.links(c) == f)[0][0] - ext_facet_domain.append(c) - ext_facet_domain.append(local_f) - ext_facet_domains = [(6, ext_facet_domain)] - - # Manually specify interior facets to integrate over - int_facet_domain = [] - for f in marked_int_facets: - if f >= facet_map.size_local or len(f_to_c.links(f)) != 2: - continue - c_0, c_1 = f_to_c.links(f)[0], f_to_c.links(f)[1] - local_f_0 = np.where(c_to_f.links(c_0) == f)[0][0] - local_f_1 = np.where(c_to_f.links(c_1) == f)[0][0] - int_facet_domain.append(c_0) - int_facet_domain.append(local_f_0) - int_facet_domain.append(c_1) - int_facet_domain.append(local_f_1) - int_facet_domains = [(3, int_facet_domain)] - - # Create measures - dx_manual = ufl.Measure("dx", subdomain_data=cell_domains, domain=msh) - ds_manual = ufl.Measure("ds", subdomain_data=ext_facet_domains, domain=msh) - dS_manual = ufl.Measure("dS", subdomain_data=int_facet_domains, domain=msh) - - # Assemble forms and check - a, L = create_forms(dx_manual, ds_manual, dS_manual) - A = fem.assemble_matrix(a) - A.scatter_reverse() - b = fem.assemble_vector(L) - - assert np.allclose(A.data, A_mt.data) - assert np.allclose(b.array, b_mt.array) \ No newline at end of file diff --git a/tests/temperer3D_mapA_example.py b/tests/temperer3D_mapA_example.py deleted file mode 100644 index 8468e73..0000000 --- a/tests/temperer3D_mapA_example.py +++ /dev/null @@ -1,214 +0,0 @@ -import numpy as np -import pickle -import itertools -from warmth3d.fixed_mesh_model import UniformNodeGridFixedSizeMeshModel -from warmth3d.Helpers import NodeGrid -from warmth3d.resqpy_helpers import read_mesh_resqml - -# -# a map for how many meters were eroded between 100-120 My -# a map for how many meters were eroded between 120-140 My => constant 50m -# - -def tic(): - #Homemade version of matlab tic and toc functions - import time - global startTime_for_tictoc - startTime_for_tictoc = time.time() - -def toc(msg=""): - import time - if 'startTime_for_tictoc' in globals(): - delta = time.time() - startTime_for_tictoc - print (msg+": Elapsed time is " + str(delta) + " seconds.") - return delta - else: - print ("Toc: start time not set") - -def nodeFromDataArray(nodeGrid, node_data_array, nodeid_x, nodeid_y): - origin_x = nodeGrid.origin_x - origin_y = nodeGrid.origin_y - stepx = nodeGrid.step_x - stepy = nodeGrid.step_y - class GenericObject(object): - pass - node = GenericObject - node.X, node.Y = node_data_array[0:2] - node.subsidence = node_data_array[3] - node.crust_ls = node_data_array[4] - node.lith_ls = node_data_array[5] - node.sed = node_data_array[2] - node.X, node.Y = (origin_x + stepx*nodeid_x, origin_y + stepy*nodeid_y) - node.sed_thickness_ls = node.sed[-1,1,:] - node.sed[0,0,:] - node.sediments = node_data_array[6] - if len(node_data_array)>7: - node.beta = node_data_array[7] - if len(node_data_array)>9: - node.depth_out = node_data_array[8].copy() - node.temperature_out = node_data_array[9].copy() - if len(node_data_array)>10: - node.parameters = node_data_array[10] - return node - -def interpolateNode(interpolationNodes, interpolationWeights): - assert len(interpolationNodes)>0 - assert len(interpolationNodes)==len(interpolationWeights) - wsum = np.sum(np.array(interpolationWeights)) - iWeightNorm = [ w/wsum for w in interpolationWeights] - class GenericObject(object): - pass - node = GenericObject - node.X = np.sum( np.array( [node.X * w for node,w in zip(interpolationNodes,iWeightNorm)] ) ) - node.Y = np.sum( np.array( [node.Y * w for node,w in zip(interpolationNodes,iWeightNorm)] ) ) - node.subsidence = np.sum( np.array( [node.subsidence * w for node,w in zip(interpolationNodes,iWeightNorm)] ) , axis = 0) - node.crust_ls = np.sum( np.array( [node.crust_ls * w for node,w in zip(interpolationNodes,iWeightNorm)] ) , axis = 0) - node.lith_ls = np.sum( np.array( [node.lith_ls * w for node,w in zip(interpolationNodes,iWeightNorm)] ) , axis = 0) - node.beta = np.sum( np.array( [node.beta * w for node,w in zip(interpolationNodes,iWeightNorm)] ) , axis = 0) - node.kAsth = np.sum( np.array( [node.kAsth * w for node,w in zip(interpolationNodes,iWeightNorm)] ) , axis = 0) - node.kLith = np.sum( np.array( [node.kLith * w for node,w in zip(interpolationNodes,iWeightNorm)] ) , axis = 0) - node.sediments = interpolationNodes[0].sediments.copy() - return node - -def loadOrInterpolateNode(nodeGrid, nodeDirectoryPrefix, nodeid_x, nodeid_y): - if (nodeid_x % 10 > 0) or (nodeid_y % 10 > 0): - nidx0 = nodeid_x - (nodeid_x % 10) - nidy0 = nodeid_y - (nodeid_y % 10) - interpolationNodes = [] - interpolationWeight = [] - sedFileName = nodeDirectoryPrefix+"node-sed-"+str(nodeid_x)+"-"+str(nodeid_y)+".pkl" - try: - sed_data_array = pickle.load( open(sedFileName,"rb") ) - X, Y = sed_data_array[0], sed_data_array[1] - sed = sed_data_array[2] - except FileNotFoundError as e: - print("1D sediment solution not found: ", sedFileName) - breakpoint() - return - ddx = [0,10] if (nodeid_x % 10 > 0) else [0] - ddy = [0,10] if (nodeid_y % 10 > 0) else [0] - for dx,dy in list(itertools.product(ddx, ddy)): - nxpos, nypos = nidx0+dx, nidy0+dy - nodeFileName = nodeDirectoryPrefix+"node-slim-"+str(nxpos)+"-"+str(nypos)+".pkl" - try: - node_data_array = pickle.load( open(nodeFileName,"rb") ) - interpolationNodes.append(nodeFromDataArray(nodeGrid,node_data_array,nxpos,nypos)) - dist = np.sqrt((nodeid_x-nxpos)**2 + (nodeid_y-nypos)**2) - interpolationWeight.append(1/dist) - except FileNotFoundError as e: - print("1D Node solution not found: ", nodeFileName) - breakpoint() - return - node = interpolateNode(interpolationNodes, interpolationWeight) - node.sed = sed - node.sed_thickness_ls = node.sed[-1,1,:] - node.sed[0,0,:] - origin_x = nodeGrid.origin_x - origin_y = nodeGrid.origin_y - stepx = nodeGrid.step_x - stepy = nodeGrid.step_y - node.X, node.Y = (origin_x + stepx*nodeid_x, origin_y + stepy*nodeid_y) - return node - else: - nodeFileName = nodeDirectoryPrefix+"node-slim-"+str(nodeid_x)+"-"+str(nodeid_y)+".pkl" - try: - node_data_array = pickle.load( open(nodeFileName,"rb") ) - except FileNotFoundError as e: - print("1D Node solution not found: ", nodeFileName) - breakpoint() - return - node = nodeFromDataArray(nodeGrid, node_data_array, nodeid_x, nodeid_y) - return node - -def run( nodeGrid, run_simulation=True, start_time=182, end_time=0, out_dir = "out-mapA/"): - nodes = [] - modelNamePrefix = ng.modelNamePrefix - - cols = ng.num_nodes_x - rows = ng.num_nodes_y - ind_step = ng.step_index - start_index_x = ng.start_index_x - start_index_y = ng.start_index_y - - convexHullEdges = [] - - for j in range(rows): - for i in range(cols): - nodeid_x = start_index_x + i * ind_step - nodeid_y = start_index_y + j * ind_step - - node = loadOrInterpolateNode(nodeGrid , nodeGrid.nodeDirectoryPrefix, nodeid_x, nodeid_y) - nodes.append(node) - - characteristic_length = ind_step * ng.step_x - - nums = 4 - dt = 314712e8 / nums - - mms2 = [] - mms_tti = [] - - tti = 0 - subvolumes = [] - - writeout = True - - if not run_simulation: - return - time_solve = 0.0 - - for tti in range(start_time, end_time-1,-1): - rebuild_mesh = (tti==start_time) - if rebuild_mesh: - print("Rebuild/reload mesh at tti=", tti) - mm2 = UniformNodeGridFixedSizeMeshModel(nodeGrid, nodes, modelName=modelNamePrefix+str(tti)) - mm2.buildMesh(tti) - else: - print("Re-generating mesh vertices at tti=", tti) - mm2.updateMesh(tti) - - print("===",tti,"=========== ") - if ( len(mms2) == 0): - tic() - mm2.setupSolverAndSolve(no_steps=40, time_step = 314712e8 * 2e2, skip_setup=False) - time_solve = time_solve + toc(msg="setup solver and solve") - else: - tic() - # mm2.setupSolverAndSolve( initial_state_model=mms2[-1], no_steps=nums, time_step=dt, skip_setup=(not rebuild_mesh)) - mm2.setupSolverAndSolve( no_steps=nums, time_step=dt, skip_setup=(not rebuild_mesh)) - time_solve = time_solve + toc(msg="setup solver and solve") - # subvolumes.append(mm2.evaluateVolumes()) - if (writeout): - tic() - mm2.writeLayerIDFunction(out_dir+"LayerID-"+str(tti)+".xdmf", tti=tti) - mm2.writeTemperatureFunction(out_dir+"Temperature-"+str(tti)+".xdmf", tti=tti) - # mm2.writeOutputFunctions(out_dir+"test4-"+str(tti)+".xdmf", tti=tti) - toc(msg="write function") - mms2.append(mm2) - mms_tti.append(tti) - print("total time solve: " , time_solve) - EPCfilename = mm2.write_hexa_mesh_resqml("temp/") - print("RESQML model written to: " , EPCfilename) - read_mesh_resqml(EPCfilename) # test reading of the .epc file - - - -# -# NOTE: to compute the required 1D node solutions, you must first run warmth3d/parallel-1Dsed.py using the same NodeGrid parameters as below! -# -ng = NodeGrid(150, 0, 485, 548, 500, 1000, 5, 1000, 1000) -ng = NodeGrid(0, 0, 11, 11, 500, 1000, 5, 100, 100) -# ng = NodeGrid(25400, 41600, 31, 31, 100, 100, 2, 100, 100) - -ng.modelNamePrefix = "mapB-121-nodes-" -ng.nodeDirectoryPrefix = "nodes-mapB/" - -import os -for output_dir in ["out-mapB/", "mesh", "temp"]: - if not os.path.exists(output_dir): - os.makedirs(output_dir) - -# import cProfile -# cProfile.run('run( ng, run_simulation=True, start_time=182, end_time=0, out_dir = "out-mapA/")') - -run( ng, run_simulation=True, start_time=170, end_time=0, out_dir = "out-mapB/") - - diff --git a/tests/ttest_from_grid.py b/tests/test_from_grid.py similarity index 93% rename from tests/ttest_from_grid.py rename to tests/test_from_grid.py index bc82bab..17209a1 100644 --- a/tests/ttest_from_grid.py +++ b/tests/test_from_grid.py @@ -1,13 +1,10 @@ -import os -import multiprocessing from pathlib import Path -import sys import warmth -maps_dir = Path("./docs/data/mapA") +maps_dir = Path("./docs/notebooks/data/") model = warmth.Model() @@ -52,7 +49,6 @@ model.simulator.run(save=False,purge=True) -# %% for i in model.builder.iter_node(): if i is not False: print(i.result.heatflow(0)) diff --git a/tests/test_grid.py b/tests/test_grid.py deleted file mode 100644 index f3188f9..0000000 --- a/tests/test_grid.py +++ /dev/null @@ -1,3 +0,0 @@ -from warmth.build import Grid -import numpy as np - diff --git a/tests/warmth3d/test_3d.py b/tests/warmth3d/test_3d.py new file mode 100644 index 0000000..1f85e38 --- /dev/null +++ b/tests/warmth3d/test_3d.py @@ -0,0 +1,97 @@ +import warmth +from pathlib import Path +from warmth.mesh_model import run_3d +import os +import numpy as np + +def test_3d_compare(): + maps_dir = Path("./docs/notebooks/data/") + model = warmth.Model() + + inputs = model.builder.input_horizons_template + + #Add surface grids to the table. You can use other method as well + inputs.loc[0]=[0,"0.gri",None,"Onlap"] + inputs.loc[1]=[66,"66.gri",None,"Onlap"] + inputs.loc[2]=[100,"100.gri",None,"Onlap"] + inputs.loc[3]=[163,"163.gri",None,"Erosive"] + inputs.loc[6]=[182,"182.gri",None,"Erosive"] + model.builder.input_horizons=inputs + + + inc = 2000 + model.builder.define_geometry(maps_dir/"0.gri",xinc=inc,yinc=inc,fformat="irap_binary") + + model.builder.extract_nodes(4,maps_dir) + + from warmth.data import haq87 + model.builder.set_eustatic_sea_level(haq87) + + for i in model.builder.iter_node(): + i.rift=[[182,175]] + # + # set 1D node parameters to be most similar to those in the (later) 3D simulation, for better comparison + # + i.bflux = False + i.adiab = 0.3e-3 + + + model.simulator.simulate_every = 1 + + # + # set 1D simulation parameters to be most similar to those in the (later) 3D simulation, for better comparison + # + model.parameters.HPdcr = 1e32 # "infinite" decay of crustal HP + model.parameters.bflux = False + model.parameters.tetha = 0 + model.parameters.alphav = 0 + + model.simulator.run(save=True,purge=True) + + + + + + try: + os.mkdir('mesh') + except FileExistsError: + pass + try: + os.mkdir('temp') + except FileExistsError: + pass + mm2 = run_3d(model.builder,model.parameters,start_time=model.parameters.time_start,end_time=0,writeout=False) + + + hx = model.builder.grid.num_nodes_x // 2 + hy = model.builder.grid.num_nodes_y // 2 + # hx = model.builder.grid.num_nodes_x - 1 - pad + # hy = model.builder.grid.num_nodes_y - 1 - pad + + nn = model.builder.nodes[hy][hx] + dd = nn._depth_out[:,0] + + node_ind = hy*model.builder.grid.num_nodes_x + hx + v_per_n = int(mm2.mesh_vertices.shape[0]/(model.builder.grid.num_nodes_y*model.builder.grid.num_nodes_x)) + + temp_1d = np.nan_to_num(nn.temperature_out[:,0], nan=5.0) + temp_3d_ind = np.array([ np.where([mm2.mesh_reindex==i])[1][0] for i in range(node_ind*v_per_n, (node_ind+1)*v_per_n) ] ) + dd_mesh = mm2.mesh.geometry.x[temp_3d_ind,2] + temp_3d_mesh = mm2.u_n.x.array[temp_3d_ind] + + temp_1d_at_mesh_pos = np.interp(dd_mesh, dd, temp_1d) + dd_subset = np.where(dd_mesh<5000) + + max_abs_error = np.amax(np.abs(temp_1d_at_mesh_pos-temp_3d_mesh)) + max_abs_error_shallow = np.amax(np.abs(temp_1d_at_mesh_pos[dd_subset]-temp_3d_mesh[dd_subset])) + print(f'Max. absolute error in temperature at 3D mesh vertex positions: {max_abs_error}') + print(f'Max. absolute error at depths < 5000m: {max_abs_error_shallow}') + + assert (max_abs_error<25.0), "Temperature difference between 3D and 1D simulations is >25" + assert (max_abs_error_shallow<5.0), "Temperature difference between 3D and 1D simulations is >5 in the sediments." + + + + +if __name__ == "__main__": + test_3d_compare() \ No newline at end of file diff --git a/tests/warmth3d/test_from_grid.py b/tests/warmth3d/test_from_grid.py deleted file mode 100644 index d16f74c..0000000 --- a/tests/warmth3d/test_from_grid.py +++ /dev/null @@ -1,135 +0,0 @@ -import warmth -from pathlib import Path - -maps_dir = Path("../../docs/notebooks/data/") -model = warmth.Model() - -inputs = model.builder.input_horizons_template - -#Add surface grids to the table. You can use other method as well -inputs.loc[0]=[0,"0.gri",None,"Onlap"] -inputs.loc[1]=[66,"66.gri",None,"Onlap"] -inputs.loc[2]=[100,"100.gri",None,"Onlap"] -inputs.loc[3]=[163,"163.gri",None,"Erosive"] -inputs.loc[4]=[168,"168.gri",None,"Erosive"] -inputs.loc[5]=[170,"170.gri",None,"Onlap"] -inputs.loc[6]=[182,"182.gri",None,"Erosive"] -model.builder.input_horizons=inputs - - -inc = 2000 -model.builder.define_geometry(maps_dir/"0.gri",xinc=inc,yinc=inc,fformat="irap_binary") - -model.builder.extract_nodes(4,maps_dir) - -from warmth.data import haq87 -model.builder.set_eustatic_sea_level(haq87) - -for i in model.builder.iter_node(): - i.rift=[[182,175]] - # - # set 1D node parameters to be most similar to those in the (later) 3D simulation, for better comparison - # - i.bflux = False - i.adiab = 0.3e-3 - # i.crustRHP = 0.0 - # i.sediments['rhp']=[0,0,0,0,0,0] - -model.simulator.simulate_every = 1 - -# -# set 1D simulation parameters to be most similar to those in the (later) 3D simulation, for better comparison -# -model.parameters.HPdcr = 1e32 # "infinite" decay of crustal HP -model.parameters.bflux = False -model.parameters.tetha = 0 -model.parameters.alphav = 0 - -model.simulator.run(save=True,purge=True) - - - -# interpolate and extrapolate the missing nodes -# find nearby nodes from the array indexer_full_sim, which is sorted by x-index -import itertools -from warmth3d.fixed_mesh_model import interpolateNode -for ni in range(len(model.builder.nodes)): - for nj in range(len(model.builder.nodes[ni])): - if model.builder.nodes[ni][nj] is False: - closest_x_up = [] - for j in range(ni,len(model.builder.nodes[nj])): - matching_x = [ i[0] for i in model.builder.indexer_full_sim if i[0]==j ] - closest_x_up = closest_x_up + list(set(matching_x)) - if len(matching_x)>0: - break - closest_x_down = [] - for j in range(ni-1,-1,-1): - matching_x = [ i[0] for i in model.builder.indexer_full_sim if i[0]==j ] - closest_x_down = closest_x_down + list(set(matching_x)) - if len(matching_x)>0: - break - closest_y_up = [] - for j in range(nj,len(model.builder.nodes[ni])): - matching_y = [ i[1] for i in model.builder.indexer_full_sim if (i[1]==j and ((i[0] in closest_x_up) or i[0] in closest_x_down)) ] - closest_y_up = closest_y_up + list(set(matching_y)) - if len(matching_y)>0: - break - closest_y_down = [] - for j in range(nj-1,-1,-1): - matching_y = [ i[1] for i in model.builder.indexer_full_sim if (i[1]==j and (i[0] in closest_x_up or i[0] in closest_x_down) ) ] - closest_y_down = closest_y_down + list(set(matching_y)) - if len(matching_y)>0: - break - - interpolationNodes = [ model.builder.nodes[i[0]][i[1]] for i in itertools.product(closest_x_up+closest_x_down, closest_y_up+closest_y_down) ] - interpolationNodes = [nn for nn in interpolationNodes if nn is not False] - node = interpolateNode(interpolationNodes) - node.X, node.Y = model.builder.grid.location_grid[ni,nj,:] - model.builder.nodes[ni][nj] = node - else: - node = interpolateNode([model.builder.nodes[ni][nj]]) # "interpolate" the node from itself to make sure the same member variables exist at the end - model.builder.nodes[ni][nj] = node - -from warmth.mesh_model import run -import os -try: - os.mkdir('mesh') -except FileExistsError: - pass -try: - os.mkdir('temp') -except FileExistsError: - pass -mm2 = run(model,start_time=model.parameters.time_start,end_time=0) - -import numpy as np -hx = model.builder.grid.num_nodes_x // 2 -hy = model.builder.grid.num_nodes_y // 2 -# hx = model.builder.grid.num_nodes_x - 1 - pad -# hy = model.builder.grid.num_nodes_y - 1 - pad - -nn = model.builder.nodes[hy][hx] -dd = nn._depth_out[:,0] - -node_ind = hy*model.builder.grid.num_nodes_x + hx -v_per_n = int(mm2.mesh_vertices.shape[0]/(model.builder.grid.num_nodes_y*model.builder.grid.num_nodes_x)) - -temp_1d = np.nan_to_num(nn.temperature_out[:,0], nan=5.0) -temp_3d_ind = np.array([ np.where([mm2.mesh_reindex==i])[1][0] for i in range(node_ind*v_per_n, (node_ind+1)*v_per_n) ] ) -dd_mesh = mm2.mesh.geometry.x[temp_3d_ind,2] -temp_3d_mesh = mm2.u_n.x.array[temp_3d_ind] - -temp_1d_at_mesh_pos = np.interp(dd_mesh, dd, temp_1d) -dd_subset = np.where(dd_mesh<5000) - -max_abs_error = np.amax(np.abs(temp_1d_at_mesh_pos-temp_3d_mesh)) -max_abs_error_shallow = np.amax(np.abs(temp_1d_at_mesh_pos[dd_subset]-temp_3d_mesh[dd_subset])) -print(f'Max. absolute error in temperature at 3D mesh vertex positions: {max_abs_error}') -print(f'Max. absolute error at depths < 5000m: {max_abs_error_shallow}') - -assert (max_abs_error<25.0), "Temperature difference between 3D and 1D simulations is >25" -assert (max_abs_error_shallow<5.0), "Temperature difference between 3D and 1D simulations is >5 in the sediments." - - - - diff --git a/warmth/build.py b/warmth/build.py index 84bafb5..3ad2044 100644 --- a/warmth/build.py +++ b/warmth/build.py @@ -668,6 +668,7 @@ def extract_nodes( for future in concurrent.futures.as_completed(futures): sed = future.result() # This will also raise any exceptions sediments_all.append(sed) + self._create_nodes(sediments_all) return @@ -681,54 +682,61 @@ def _create_nodes(self, all_sediments_grid: List[List[List]]): Extracted sediment objects """ indexer = self.grid.indexing_arr + + for index in indexer: + if self.nodes[index[0]][index[1]] != False: node_sed: list[_sediment_layer_] = [] for sed_grid in all_sediments_grid: node_sed.append(sed_grid[index[0]][index[1]]) - top = np.empty(0) - topage = np.empty(0) - k_cond = np.empty(0) - rhp = np.empty(0) - phi = np.empty(0) - decay = np.empty(0) - solidus = np.empty(0) - liquidus = np.empty(0) - strat = np.empty(0,dtype=str) - inputRef = np.empty(0,dtype=int) - for hor in node_sed: - top = np.append(top, float(hor.top)) - topage = np.append(topage, int(hor.topage)) - k_cond = np.append(k_cond, float(hor.thermoconductivity)) - rhp = np.append(rhp, float(hor.rhp)) - phi = np.append(phi, float(hor.phi)) - decay = np.append(decay, float(hor.decay)) - solidus = np.append(solidus, float(hor.solidus)) - liquidus = np.append(liquidus, float(hor.liquidus)) - strat = np.append(strat,hor.strat) - inputRef= np.append(inputRef,hor.horizon_index) - - - df = pd.DataFrame({'top': top, 'topage': topage, 'k_cond': k_cond, - 'rhp': rhp, 'phi': phi, 'decay': decay, 'solidus': solidus, 'liquidus': liquidus,'strat':strat,'horizonIndex':inputRef}) - df = df.sort_values(by=["topage"],ignore_index=True) - - - #df.reset_index(drop=True,inplace=True) - df.at[2, 'top'] = np.nan - - df.at[3, 'top'] = np.nan - checker = self._check_nan_sed(df) - if checker is False: - self.nodes[index[0]][index[1]] = False + if all(node_sed) is False: + logger.warning(f"dropping node {index}. One of the layer has no dpeth value") else: - df = self._fix_nan_sed(df) - n = single_node() - n.X=node_sed[0].X - n.Y=node_sed[0].Y - n.sediments_inputs=df - n.indexer = index - self.nodes[index[0]][index[1]] = n + top = np.empty(0) + topage = np.empty(0) + k_cond = np.empty(0) + rhp = np.empty(0) + phi = np.empty(0) + decay = np.empty(0) + solidus = np.empty(0) + liquidus = np.empty(0) + strat = np.empty(0,dtype=str) + inputRef = np.empty(0,dtype=int) + for hor in node_sed: + top = np.append(top, float(hor.top)) + topage = np.append(topage, int(hor.topage)) + k_cond = np.append(k_cond, float(hor.thermoconductivity)) + rhp = np.append(rhp, float(hor.rhp)) + phi = np.append(phi, float(hor.phi)) + decay = np.append(decay, float(hor.decay)) + solidus = np.append(solidus, float(hor.solidus)) + liquidus = np.append(liquidus, float(hor.liquidus)) + strat = np.append(strat,hor.strat) + inputRef= np.append(inputRef,hor.horizon_index) + + + df = pd.DataFrame({'top': top, 'topage': topage, 'k_cond': k_cond, + 'rhp': rhp, 'phi': phi, 'decay': decay, 'solidus': solidus, 'liquidus': liquidus,'strat':strat,'horizonIndex':inputRef}) + df = df.sort_values(by=["topage"],ignore_index=True) + + + #df.reset_index(drop=True,inplace=True) + df.at[2, 'top'] = np.nan + + df.at[3, 'top'] = np.nan + checker = self._check_nan_sed(df) + + if checker is False: + self.nodes[index[0]][index[1]] = False + else: + df = self._fix_nan_sed(df) + n = single_node() + n.X=node_sed[0].X + n.Y=node_sed[0].Y + n.sediments_inputs=df + n.indexer = index + self.nodes[index[0]][index[1]] = n else: pass return diff --git a/warmth/forward_modelling.py b/warmth/forward_modelling.py index baef0a3..d461b3b 100644 --- a/warmth/forward_modelling.py +++ b/warmth/forward_modelling.py @@ -1586,11 +1586,7 @@ def calculate_new_temperature(self, mean_porosity_arr, self.current_node.sediments["solidus"].values[sed_idx_arr]) conductivity_sed = self._sediment_conductivity_sekiguchi( mean_porosity_arr, self.current_node.sediments["k_cond"].values[sed_idx_arr],Tsed) - if (self.current_node.X==12150) and (self.current_node.Y==12000): - print("conductivity_sed", conductivity_sed[0:20], conductivity_sed[-20:] ) - # print("conductivity_sed", np.mean(conductivity_sed), np.nanmin(conductivity_sed),np.nanmax(conductivity_sed) ) - # print("conductivity_sed", np.mean(xsed), np.nanmin(xsed),np.nanmax(xsed) ) - print("conductivity_sed", len(conductivity_sed), len(xsed)) + conductivity_crust_lith = self._build_crust_lithosphere_properties( coord_crust_lith, hc, hLith, self.current_node.kCrust, self.current_node.kLith, self.current_node.kAsth) diff --git a/warmth/mesh_model.py b/warmth/mesh_model.py index 4dfd636..f89f661 100644 --- a/warmth/mesh_model.py +++ b/warmth/mesh_model.py @@ -1,3 +1,4 @@ +from typing import Tuple import numpy as np from mpi4py import MPI import meshio @@ -5,11 +6,11 @@ from petsc4py import PETSc import ufl from scipy.interpolate import LinearNDInterpolator - -from warmth.build import single_node -from .model import Model +from progress.bar import Bar +from warmth.build import single_node, Builder +from .parameters import Parameters from warmth.logging import logger -from .mesh_utils import top_crust,top_sed,thick_crust, top_lith, top_asth, top_sed_id, bottom_sed_id,NodeGrid +from .mesh_utils import top_crust,top_sed,thick_crust, top_lith, top_asth, top_sed_id, bottom_sed_id,interpolate_all_nodes from .resqpy_helpers import write_tetra_grid_with_properties, write_hexa_grid_with_properties,read_mesh_resqml_hexa def tic(): #Homemade version of matlab tic and toc functions @@ -36,12 +37,14 @@ class UniformNodeGridFixedSizeMeshModel: point_domain_edge_map = {} point_top_vertex_map = {} point_bottom_vertex_map = {} - def __init__(self, model:Model, modelName="test", sedimentsOnly = False): - self.node1D = [n for n in model.builder.iter_node()] + def __init__(self, builder:Builder,parameters:Parameters, sedimentsOnly = False): + self._builder = builder + self._parameters=parameters + self.node1D = [n for n in self._builder.iter_node()] self.num_nodes = len(self.node1D) self.mesh = None - self.modelName = modelName + self.modelName = self._parameters.name self.Temp0 = 5 self.TempBase = 1369 self.verbose = True @@ -55,8 +58,8 @@ def __init__(self, model:Model, modelName="test", sedimentsOnly = False): self.numElemInAsth = 0 if self.runSedimentsOnly else 2 # split asth hexahedron into pieces - self.num_nodes_x = model.builder.grid.num_nodes_x - self.num_nodes_y = model.builder.grid.num_nodes_y + self.num_nodes_x = self._builder.grid.num_nodes_x + self.num_nodes_y = self._builder.grid.num_nodes_y self.convexHullEdges = [] for i in range(self.num_nodes_x-1): edge = [i, i+1] @@ -79,7 +82,7 @@ def __init__(self, model:Model, modelName="test", sedimentsOnly = False): self.thermalCond = None self.mean_porosity = None self.c_rho = None - self.numberOfSediments = model.builder.input_horizons.shape[0]-1 #skip basement + self.numberOfSediments = self._builder.input_horizons.shape[0]-1 #skip basement self.numberOfSedimentCells = self.numberOfSediments * self.numElemPerSediment self.interpolators = {} @@ -260,14 +263,14 @@ def getTopOfAsthAtNode(self, tti, node:single_node): return z0 # - def getSedimentPropForLayerID(self, property, layer_id, node_index): + def getSedimentPropForLayerID(self, property, layer_id:int, node_index:int) ->float: """ """ assert property in ['k_cond', 'rhp', 'phi', 'decay', 'solidus', 'liquidus'], "Unknown property " + property if (layer_id>=0) and (layer_id0.7 and node.Y<40000: # print("phi", property, phi, node_index, node) # breakpoint() @@ -276,7 +279,7 @@ def getSedimentPropForLayerID(self, property, layer_id, node_index): # breakpoint() # phi = self.globalSediments[property][layer_id] # assert abs(phi-phi0)<1e-6 - return phi + return prop if (layer_id<=-1) and (layer_id>=-3): lid = -layer_id -1 node = self.node1D[node_index] @@ -294,7 +297,7 @@ def getSedimentPropForLayerID(self, property, layer_id, node_index): return [node.crustliquid,node.lithliquid,node.asthliquid][lid] # liquid density for crust, lith, aest return np.nan - def porosity0ForLayerID(self, layer_id, node_index): + def porosity0ForLayerID(self, layer_id:int, node_index:int)->Tuple[float, float]: """Porosity (at surface) conductivity value for the given layer index """ if (layer_id==-1): @@ -311,37 +314,24 @@ def porosity0ForLayerID(self, layer_id, node_index): return phi, decay return 0.0, 0.0 - def cRhoForLayerID(self, ss, node_index): - # - # prefactor 1000 is the heat capacity.. assumed constant - # + def cRhoForLayerID(self, ss:int, node_index:int)->float: node = self.node1D[node_index] if (ss==-1): - return 1000*node.crustsolid + return self._parameters.cp*node.crustsolid if (ss==-2): - return 1000*node.lithsolid + return self._parameters.cp*node.lithsolid if (ss==-3): - return 1000*node.asthsolid + return self._parameters.cp*node.asthsolid if (ss>=0) and (ssfloat: """Thermal conductivity for a layer ID index """ - # ind0 = cell_id - # ind1 = self.cell_index[cell_id] - # ind2 = np.where(self.cell_index==cell_id)[0][0] - # if self.cell_data_layerID[ind2] != ss: - # print("ind", cell_id, ind0, ind1, ind2) - # print(len(self.cell_index), np.amin(self.cell_index), np.amax(self.cell_index)) - # print(len(self.cell_data_layerID), np.amin(self.cell_data_layerID), np.amax(self.cell_data_layerID)) - # print("kForLayerID", ss, ind0, ind1, ind2, self.cell_data_layerID[self.cell_index[cell_id]] ) - # assert self.cell_data_layerID[ind2] == ss if (node_index > len(self.node1D)-1): - print("cell ID", node_index, len(self.node1D)) - breakpoint() + raise Exception(f"Node index {node_index} larger then node length {len(self.node1D)}") node = self.node1D[node_index] if (ss==-1): return node.kCrust @@ -356,17 +346,13 @@ def kForLayerID(self, ss, node_index): return kc - def rhpForLayerID(self, ss, node_index): + def rhpForLayerID(self, ss:int, node_index:int)->float: """Radiogenic heat production for a layer ID index """ if (ss==-1): node = self.node1D[node_index] kc = node.crustRHP * node._upperCrust_ratio return kc - elif (ss==-2): - return 0 - elif (ss==-3): - return 0 elif (ss>=0) and (ss0) else 3600*24*365 * 5000000 - num_steps = no_steps + num_steps = n_steps # # solver setup, see: @@ -929,18 +908,16 @@ def setupSolverAndSolve(self, time_step=-1, no_steps=100, skip_setup = False, in v = ufl.TestFunction(self.V) a = self.c_rho*u*v*ufl.dx + dt*ufl.dot(self.thermalCond*ufl.grad(u), ufl.grad(v)) * ufl.dx + f = self.rhpFcn + logger.info(f"mean RHP {np.mean(self.rhpFcn.x.array[:])}") - # source = self.globalSediments.rhp[self.numberOfSediments-1] * 1e-6 # conversion from uW/m^3 - # f = dolfinx.fem.Constant(self.mesh, PETSc.ScalarType(source)) # source term - f = self.rhpFcn # * 1e-6 # conversion from uW/m^3 - print("mean RHP", np.mean(self.rhpFcn.x.array[:])) if ( self.useBaseFlux ): # baseFlux = 0.03 if (self.tti>50) else 0.03 baseFlux = self.baseFluxMagnitude # define Neumann condition: constant flux at base # expression g defines values of Neumann BC (heat flux at base) - x = ufl.SpatialCoordinate(self.mesh) + #x = ufl.SpatialCoordinate(self.mesh) domain_c = dolfinx.fem.Function(self.V) if (self.CGorder>1): def marker(x): @@ -948,15 +925,17 @@ def marker(x): return x[2,:]>3990 facets = dolfinx.mesh.locate_entities_boundary(self.mesh, dim=(self.mesh.topology.dim - 2), marker=marker ) - print(type(facets), facets.shape) + #print(type(facets), facets.shape) dofs = dolfinx.fem.locate_dofs_topological(V=self.V, entity_dim=1, entities=facets) - print( type(dofs), len(dofs)) - print(facets.shape, dofs.shape) + #print( type(dofs), len(dofs)) + #print(facets.shape, dofs.shape) if (len(facets)>0): - print( np.amax(facets)) + #print( np.amax(facets)) + pass if (len(dofs)>0): - print( np.amax(dofs)) - print(type(domain_c.x.array), len(domain_c.x.array)) + #print( np.amax(dofs)) + pass + #print(type(domain_c.x.array), len(domain_c.x.array)) domain_c.x.array[ dofs ] = 1 else: basepos = self.getBaseAtMultiplePos(self.mesh.geometry.x[:,0], self.mesh.geometry.x[:,1]) @@ -1298,7 +1277,7 @@ def interpolateResult(self, x): print("PING V V", point) def boundary(x): return np.full(x.shape[1], True) - entities = dolfinx.mesh.locate_entities(self.mesh, 3, boundary ) + #entities = dolfinx.mesh.locate_entities(self.mesh, 3, boundary ) breakpoint() points_on_proc.append(point) points_cells.append(cell_candidates.links(i)[0]) @@ -1330,57 +1309,49 @@ def pointIsOnDomainEdge(self, pt, node0, node1, weight): return False -def run( model:Model, run_simulation=True, start_time=182, end_time=0, out_dir = "out-mapA/"): - - +def run_3d( builder:Builder, parameters:Parameters, start_time=182, end_time=0, out_dir = "out-mapA/",sedimentsOnly=False,writeout=True): + builder=interpolate_all_nodes(builder) nums = 4 - dt = 314712e8 / nums - + dt = parameters.myr2s / nums # time step is 1/4 of 1Ma mms2 = [] mms_tti = [] - tti = 0 - subvolumes = [] - - writeout = True - if not run_simulation: - return time_solve = 0.0 - - for tti in range(start_time, end_time-1,-1): #start from oldest - rebuild_mesh = (tti==start_time) - if rebuild_mesh: - print("Rebuild/reload mesh at tti=", tti) - mm2 = UniformNodeGridFixedSizeMeshModel(model, modelName="test"+str(tti)) - print("builing") - mm2.buildMesh(tti) - print("done") - else: - print("Re-generating mesh vertices at tti=", tti) - mm2.updateMesh(tti) - - print("===",tti,"=========== ") - if ( len(mms2) == 0): - tic() - mm2.setupSolverAndSolve(no_steps=40, time_step = 314712e8 * 2e2, skip_setup=False) - time_solve = time_solve + toc(msg="setup solver and solve") - else: - tic() - # mm2.setupSolverAndSolve( initial_state_model=mms2[-1], no_steps=nums, time_step=dt, skip_setup=(not rebuild_mesh)) - mm2.setupSolverAndSolve( no_steps=nums, time_step=dt, skip_setup=(not rebuild_mesh)) - time_solve = time_solve + toc(msg="setup solver and solve") - # subvolumes.append(mm2.evaluateVolumes()) - if (writeout): - tic() - mm2.writeLayerIDFunction(out_dir+"LayerID-"+str(tti)+".xdmf", tti=tti) - mm2.writeTemperatureFunction(out_dir+"Temperature-"+str(tti)+".xdmf", tti=tti) - # mm2.writeOutputFunctions(out_dir+"test4-"+str(tti)+".xdmf", tti=tti) - toc(msg="write function") - mms2.append(mm2) - mms_tti.append(tti) + with Bar('Processing...',check_tty=False, max=(start_time-end_time)) as bar: + for tti in range(start_time, end_time-1,-1): #start from oldest + rebuild_mesh = (tti==start_time) + if rebuild_mesh: + logger.info(f"Rebuild/reload mesh at {tti}") + mm2 = UniformNodeGridFixedSizeMeshModel(builder, parameters,sedimentsOnly) + mm2.buildMesh(tti) + else: + logger.info(f"Re-generating mesh vertices at {tti}") + mm2.updateMesh(tti) + logger.info(f"Solving {tti}") + + if ( len(mms2) == 0): + tic() + mm2.setupSolverAndSolve(n_steps=40, time_step = 314712e8 * 2e2, skip_setup=False) + time_solve = time_solve + toc(msg="setup solver and solve") + else: + tic() + mm2.setupSolverAndSolve( n_steps=nums, time_step=dt, skip_setup=(not rebuild_mesh)) + time_solve = time_solve + toc(msg="setup solver and solve") + if (writeout): + tic() + mm2.writeLayerIDFunction(out_dir+"LayerID-"+str(tti)+".xdmf", tti=tti) + mm2.writeTemperatureFunction(out_dir+"Temperature-"+str(tti)+".xdmf", tti=tti) + # mm2.writeOutputFunctions(out_dir+"test4-"+str(tti)+".xdmf", tti=tti) + toc(msg="write function") + + mms2.append(mm2) + mms_tti.append(tti) + logger.info(f"Simulated time step {tti}") + bar.next() print("total time solve: " , time_solve) - EPCfilename = mm2.write_hexa_mesh_resqml("temp/") - print("RESQML model written to: " , EPCfilename) - read_mesh_resqml_hexa(EPCfilename) # test reading of the .epc file + if (writeout): + EPCfilename = mm2.write_hexa_mesh_resqml("temp/") + print("RESQML model written to: " , EPCfilename) + read_mesh_resqml_hexa(EPCfilename) # test reading of the .epc file return mm2 \ No newline at end of file diff --git a/warmth/mesh_utils.py b/warmth/mesh_utils.py index 53b8159..ded61f8 100644 --- a/warmth/mesh_utils.py +++ b/warmth/mesh_utils.py @@ -1,6 +1,8 @@ from dataclasses import dataclass - -from warmth.build import single_node +from typing import List +import itertools +import numpy as np +from .build import Builder, single_node @dataclass class NodeGrid: @@ -110,3 +112,81 @@ def volumeOfTet(points): cd = points[2]-points[3] bdcd = np.cross(bd,cd) return np.linalg.norm(np.dot(ad,bdcd))/6 + +def interpolateNode(interpolationNodes: List[single_node], interpolationWeights=None) -> single_node: + assert len(interpolationNodes)>0 + if interpolationWeights is None: + interpolationWeights = np.ones([len(interpolationNodes),1]) + assert len(interpolationNodes)==len(interpolationWeights) + wsum = np.sum(np.array(interpolationWeights)) + iWeightNorm = [ w/wsum for w in interpolationWeights] + + node = single_node() + node.__dict__.update(interpolationNodes[0].__dict__) + node.X = np.sum( np.array( [node.X * w for node,w in zip(interpolationNodes,iWeightNorm)] ) ) + node.Y = np.sum( np.array( [node.Y * w for node,w in zip(interpolationNodes,iWeightNorm)] ) ) + + times = range(node.result._depth.shape[1]) + if node.subsidence is None: + node.subsidence = np.sum( np.array( [ [node.result.seabed(t) for t in times] * w for node,w in zip(interpolationNodes,iWeightNorm)] ) , axis = 0) + if node.crust_ls is None: + node.crust_ls = np.sum( np.array( [ [node.result.crust_thickness(t) for t in times] * w for node,w in zip(interpolationNodes,iWeightNorm)] ) , axis = 0) + if node.lith_ls is None: + node.lith_ls = np.sum( np.array( [ [node.result.lithosphere_thickness(t) for t in times] * w for node,w in zip(interpolationNodes,iWeightNorm)] ) , axis = 0) + + if node.beta is None: + node.beta = np.sum( np.array( [node.beta * w for node,w in zip(interpolationNodes,iWeightNorm)] ) , axis = 0) + if node.kAsth is None: + node.kAsth = np.sum( np.array( [node.kAsth * w for node,w in zip(interpolationNodes,iWeightNorm)] ) , axis = 0) + if node.kLith is None: + node.kLith = np.sum( np.array( [node.kLith * w for node,w in zip(interpolationNodes,iWeightNorm)] ) , axis = 0) + if node._depth_out is None: + node._depth_out = np.sum([node.result._depth_out*w for n,w in zip(interpolationNodes[0:1], [1] )], axis=0) + if node.temperature_out is None: + node.temperature_out = np.sum([n.result.temperature_out*w for n,w in zip(interpolationNodes[0:1], [1] )], axis=0) + + if node.sed is None: + node.sed = np.sum([n.sed*w for n,w in zip(interpolationNodes,iWeightNorm)], axis=0) + if node.sed_thickness_ls is None: + node.sed_thickness_ls = node.sed[-1,1,:] - node.sed[0,0,:] + return node + + +def interpolate_all_nodes(builder:Builder)->Builder: + for ni in range(len(builder.nodes)): + for nj in range(len(builder.nodes[ni])): + if builder.nodes[ni][nj] is False: + closest_x_up = [] + for j in range(ni,len(builder.nodes[nj])): + matching_x = [ i[0] for i in builder.indexer_full_sim if i[0]==j ] + closest_x_up = closest_x_up + list(set(matching_x)) + if len(matching_x)>0: + break + closest_x_down = [] + for j in range(ni-1,-1,-1): + matching_x = [ i[0] for i in builder.indexer_full_sim if i[0]==j ] + closest_x_down = closest_x_down + list(set(matching_x)) + if len(matching_x)>0: + break + closest_y_up = [] + for j in range(nj,len(builder.nodes[ni])): + matching_y = [ i[1] for i in builder.indexer_full_sim if (i[1]==j and ((i[0] in closest_x_up) or i[0] in closest_x_down)) ] + closest_y_up = closest_y_up + list(set(matching_y)) + if len(matching_y)>0: + break + closest_y_down = [] + for j in range(nj-1,-1,-1): + matching_y = [ i[1] for i in builder.indexer_full_sim if (i[1]==j and (i[0] in closest_x_up or i[0] in closest_x_down) ) ] + closest_y_down = closest_y_down + list(set(matching_y)) + if len(matching_y)>0: + break + + interpolationNodes = [ builder.nodes[i[0]][i[1]] for i in itertools.product(closest_x_up+closest_x_down, closest_y_up+closest_y_down) ] + interpolationNodes = [nn for nn in interpolationNodes if nn is not False] + node = interpolateNode(interpolationNodes) + node.X, node.Y = builder.grid.location_grid[ni,nj,:] + builder.nodes[ni][nj] = node + else: + node = interpolateNode([builder.nodes[ni][nj]]) # "interpolate" the node from itself to make sure the same member variables exist at the end + builder.nodes[ni][nj] = node + return builder \ No newline at end of file diff --git a/warmth/parameters.py b/warmth/parameters.py index b9dd4f7..5813750 100644 --- a/warmth/parameters.py +++ b/warmth/parameters.py @@ -45,7 +45,7 @@ def __init__(self) -> None: self.maxContLith: float = 130000.0 self.starting_beta: float = 1.1 self.positive_down = True - + self.name:str = "model" pass @property diff --git a/warmth/simulator.py b/warmth/simulator.py index 65c61f7..0eee92e 100644 --- a/warmth/simulator.py +++ b/warmth/simulator.py @@ -1,4 +1,4 @@ -from multiprocessing import Pool, get_context +from multiprocessing import get_context import concurrent.futures from progress.bar import Bar from pathlib import Path