diff --git a/.github/workflows/build_documentation.yml b/.github/workflows/build_documentation.yml index 1fc28dbf1..021482be7 100644 --- a/.github/workflows/build_documentation.yml +++ b/.github/workflows/build_documentation.yml @@ -51,6 +51,7 @@ jobs: - name: Sphinx build env: PYTHON_PACKAGE_MANAGER: ${{steps.environment.outputs.ppm}} + PYTHONPATH: ${{ github.workspace }} run: | set -euxo pipefail $PYTHON_PACKAGE_MANAGER activate queens diff --git a/tests/input_files/third_party/fourc/coarse_plate_dirichlet_template.4C.yaml b/tests/input_files/third_party/fourc/coarse_plate_dirichlet_template.4C.yaml index 02e152be6..e4a568a75 100644 --- a/tests/input_files/third_party/fourc/coarse_plate_dirichlet_template.4C.yaml +++ b/tests/input_files/third_party/fourc/coarse_plate_dirichlet_template.4C.yaml @@ -5,11 +5,13 @@ PROBLEM TYPE: PROBLEMTYPE: "Structure" IO: OUTPUT_SPRING: true + STRUCT_STRESS: Cauchy IO/RUNTIME VTK OUTPUT: INTERVAL_STEPS: 1 IO/RUNTIME VTK OUTPUT/STRUCTURE: OUTPUT_STRUCTURE: true DISPLACEMENT: true + STRESS_STRAIN: true STRUCTURAL DYNAMIC: INT_STRATEGY: "Standard" RESEVERYERGY: 1 diff --git a/tutorials/3_orchestrating_4c_simulations/3_orchestrating_4c_simulations.ipynb b/tutorials/3_orchestrating_4c_simulations/3_orchestrating_4c_simulations.ipynb index 3ae6dd058..33fb940e7 100644 --- a/tutorials/3_orchestrating_4c_simulations/3_orchestrating_4c_simulations.ipynb +++ b/tutorials/3_orchestrating_4c_simulations/3_orchestrating_4c_simulations.ipynb @@ -60,7 +60,7 @@ "source": [ "# Define some paths\n", "from pathlib import Path\n", - "from utils import find_repo_root\n", + "from tutorials.utils import find_repo_root\n", "\n", "home = Path.home()\n", "NOTEBOOK_DIR = Path.cwd()\n", @@ -277,7 +277,7 @@ "source": [ "# This is just needed for some nice plots\n", "\n", - "from utils import plot_results\n", + "from plot_results import plot_results\n", "if pyvista_available:\n", " plotter = pv.Plotter()\n", " plot_results(\n", diff --git a/tutorials/3_orchestrating_4c_simulations/experimental_data.py b/tutorials/3_orchestrating_4c_simulations/experimental_data.py index 7bd67c438..62d0cdc9d 100644 --- a/tutorials/3_orchestrating_4c_simulations/experimental_data.py +++ b/tutorials/3_orchestrating_4c_simulations/experimental_data.py @@ -12,7 +12,7 @@ # should have received a copy of the GNU Lesser General Public License along with QUEENS. If not, # see . # -"""Data for tutorial 4.""" +"""Data for tutorial 3.""" import numpy as np experimental_elementwise_cauchy_stress = np.array( diff --git a/tutorials/3_orchestrating_4c_simulations/utils.py b/tutorials/3_orchestrating_4c_simulations/plot_results.py similarity index 79% rename from tutorials/3_orchestrating_4c_simulations/utils.py rename to tutorials/3_orchestrating_4c_simulations/plot_results.py index 6ddeb74aa..11a720716 100644 --- a/tutorials/3_orchestrating_4c_simulations/utils.py +++ b/tutorials/3_orchestrating_4c_simulations/plot_results.py @@ -12,9 +12,8 @@ # should have received a copy of the GNU Lesser General Public License along with QUEENS. If not, # see . # -"""Utility functions for tutorial 4.""" +"""Plotting functions for tutorial 3.""" import os -from pathlib import Path os.environ["PYVISTA_OFF_SCREEN"] = "true" os.environ["VTK_DEFAULT_RENDER_WINDOW_OFFSCREEN"] = "1" @@ -55,18 +54,3 @@ def plot_results( (0.0, 0.0, 0.2749999999999999), (-0.97637930372977, -0.08995062285804697, 0.19644933366041056), ] - - -def find_repo_root(start: Path) -> Path: - """Find the repository root by going upwards from a starting directory. - - Args: - start (Path): Starting directory to begin the upward search. - - Returns: - Path: Path to the detected repository root (containing ``pyproject.toml``). - """ - for p in [start, *start.parents]: - if (p / "pyproject.toml").exists(): - return p - raise RuntimeError(f"Could not find repo root from {start}") diff --git a/tutorials/4_quantifying_uncertainty_due_to_heterogeneous_material_fields/4_quantifying_uncertainty_due_to_heterogeneous_material_fields.ipynb b/tutorials/4_quantifying_uncertainty_due_to_heterogeneous_material_fields/4_quantifying_uncertainty_due_to_heterogeneous_material_fields.ipynb new file mode 100644 index 000000000..18e5cea93 --- /dev/null +++ b/tutorials/4_quantifying_uncertainty_due_to_heterogeneous_material_fields/4_quantifying_uncertainty_due_to_heterogeneous_material_fields.ipynb @@ -0,0 +1,616 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Quantifying uncertainty due to heterogeneous material fields\n", + "\n", + "In this tutorial we'll have a look at prescribing heterogenous material parameter fields using [4C](https://github.com/4C-multiphysics/4C) in order to do an uncertainty quantification on the resulting Cauchy stresses." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## The 4C model\n", + "\n", + "In this example, we'll have a look at a nonlinear solid mechanics problem. The momentum equation is given by\n", + "\n", + "$$\\ddot{\\boldsymbol{u}} = \\nabla \\cdot \\boldsymbol{\\sigma} + \\boldsymbol{b} \\text{ in } \\Omega \\times T$$\n", + "\n", + "where $u$ are the displacements, $\\boldsymbol{\\sigma}$ the Cauchy stress and $b$ volumetric body forces. For the examples, we assume $\\boldsymbol{b}=\\boldsymbol{0}$, and we prescribe Dirichlet boundary conditions\n", + "$$ u_x = 2.5 \\frac{t}{T_s} \\text{ on } \\Gamma_r = \\{x=12.5\\} \\times T$$\n", + "$$ u_y = 2.5 \\frac{t}{T_s} \\text{ on } \\Gamma_t = \\{y=12.5\\} \\times T$$\n", + "$$ u_y = 0 \\text{ on } \\Gamma_b = \\{y=-12.5\\} \\times T$$\n", + "$$ u_x = 0 \\text{ on } \\Gamma_l = \\{x=-12.5\\} \\times T$$\n", + "\n", + "where $T_s = 40$ and $T=[0, T_s]$. For $t=0$ we set $\\boldsymbol{u}=\\boldsymbol{0}$.\n", + "\n", + "The domain $\\Omega$ is a two-dimensional hyperelastic isochoric membrane, assuming plane stress with strain energy function\n", + "\n", + "$$\\Psi = \\frac{\\mu}{2} (J^{-2/3}I_{\\boldsymbol{C}}-3)$$\n", + "\n", + "where $\\mu$ is the shear modulus, $I_{\\boldsymbol{C}}$ is the first invariant of the Cauchy-Green tensor and $J$ the determinant of the deformation gradient. \n", + "\n", + "The problem is discretized in space using linear hexahedral finite elements and finite differences in time.\n", + "\n", + "> Note: This tutorial is based on [this](https://github.com/queens-py/queens/blob/main/tests/input_files/third_party/fourc/coarse_plate_dirichlet_template.4C.yaml) QUEENS test." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Uncertainty quantification\n", + "\n", + "The goal of this example is to compute the expectation of the Cauchy stress tensor for $t=T_s$\n", + "\n", + "$$\\overline{\\boldsymbol{\\sigma}}=\\int \\boldsymbol{\\sigma}(\\mu)p(\\mu)d\\mu$$\n", + "\n", + "where the shear modulus $\\mu$ is a random field. This scenario lends itself to study the effect of the heterogeneous constitutive material field on the Cauchy stress.\n", + "\n", + "Using the law of the unconscious statistician, we can reformulate the problem as\n", + "$$\\overline{\\boldsymbol{\\sigma}}=\\int \\boldsymbol{\\sigma}(\\mu)p(\\mu)d\\mu=\\int \\boldsymbol{\\sigma}(\\mu(\\boldsymbol{\\theta}))p(\\boldsymbol{\\theta})d\\boldsymbol{\\theta}$$\n", + "\n", + "where $\\theta$ are parameters describing the heterogeneous material field $\\mu(\\theta)$." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## The random field\n", + "\n", + "Let's define a random field:\n", + "$$\\mu(\\theta) = 0.1+\\exp\\left(-0.5\\frac{(x-\\theta_1)^2 + (y-\\theta_2)^2}{16}\\right)$$\n", + "\n", + "We set the distribution of the latent variable to be a multivariate normal distribution $p(\\boldsymbol{\\theta})=\\mathcal{N}(\\boldsymbol{\\theta}|\\boldsymbol{0}, 5\\boldsymbol{I})$ with mean value $\\boldsymbol{0}$ and a variance of $5$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "\n", + "def random_field_function(theta, positions):\n", + " diff = positions.copy()[:, :2]\n", + "\n", + " diff[:, 0] = diff[:, 0] - theta[0]\n", + " diff[:, 1] = diff[:, 1] - theta[1]\n", + " return 0.1 + np.exp(-0.5 * (np.sum(diff**2, axis=1)) / 16)\n", + "\n", + "\n", + "from queens.distributions import Normal\n", + "\n", + "theta_distribution = Normal(np.array([0, 0]), 5 * np.eye(2))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Monte Carlo integration\n", + "\n", + "Since the integral above cannot be computed analytically, we employ Monte Carlo integration\n", + "\n", + "$$\\overline{\\boldsymbol{\\sigma}}=\\int \\boldsymbol{\\sigma}(\\mu(\\boldsymbol{\\theta}))p(\\boldsymbol{\\theta})d\\boldsymbol{\\theta} \\approx \\frac{1}{N} \\sum_{s=0}^N \\boldsymbol{\\sigma}(\\mu(\\boldsymbol{\\theta}^{(s)}))$$\n", + "\n", + "where $\\boldsymbol{\\theta}^{(s)}$ are independent and identically distributed samples drawn from $p(\\boldsymbol{\\theta})$. For this examples we'll set $N=100$, so 100 4C runs will be done per QUEENS run." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Define some paths\n", + "from pathlib import Path\n", + "\n", + "try:\n", + " from tutorials.utils import find_repo_root\n", + "except ModuleNotFoundError:\n", + " # Quick fix to allow the import in a local setting\n", + " import sys\n", + " REPO_ROOT = Path.cwd().parents[1]\n", + " sys.path.insert(0, str(REPO_ROOT))\n", + " from tutorials.utils import find_repo_root\n", + "\n", + "home = Path.home()\n", + "NOTEBOOK_DIR = Path.cwd()\n", + "\n", + "QUEENS_BASE_DIR = find_repo_root(NOTEBOOK_DIR)\n", + "QUEENS_EXPERIMENTS_DIR = home / \"queens-experiments\"\n", + "TUTORIAL_DIR = QUEENS_BASE_DIR / \"tutorials\" / \"4_quantifying_uncertainty_due_to_heterogeneous_material_fields\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Let's get the domain from the 4C input file\n", + "from queens_interfaces.fourc.random_material_preprocessor import (\n", + " create_jinja_json_template,\n", + ")\n", + "import pyvista as pv\n", + "\n", + "# Open the mesh for the random field\n", + "mesh = pv.read(TUTORIAL_DIR / \"membrane_20.e\")[0][0]\n", + "\n", + "# Compute the cell centers\n", + "centers = mesh.cell_centers().points\n", + "\n", + "# Construct the RF information\n", + "mu_rf_parameters = {\n", + " \"coords\": centers,\n", + " \"keys\": [f\"MUE_{i}\" for i in range(1, 1 + len(centers))],\n", + "}\n", + "\n", + "# Contains the random field values\n", + "material_file_template = \"material.json\"\n", + "\n", + "# Create the material file for 4C\n", + "create_jinja_json_template(\n", + " \"MUE\",\n", + " np.arange(1, len(centers) + 1),\n", + " mu_rf_parameters[\"keys\"],\n", + " material_file_template,\n", + ")\n", + "\n", + "\n", + "from random_field import CustomRandomField\n", + "\n", + "random_field = CustomRandomField(\n", + " mu_rf_parameters, # Parameter names and coordinates\n", + " theta_distribution, # Latent variables\n", + " random_field_function, # Expansion, transformation from theta to mu\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now that we defined a random field, let's look at some samples!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Let's plot some samples\n", + "\n", + "import matplotlib.pyplot as plt\n", + "from plot_input_random_field import plot_field\n", + "\n", + "latent_samples = random_field.draw(3)\n", + "mu_samples = random_field.expanded_representation(latent_samples)\n", + "\n", + "fig = plt.figure(figsize=(12, 4))\n", + "\n", + "for i, mu_s in enumerate(mu_samples):\n", + " ax = fig.add_subplot(1, 3, i + 1, projection=\"3d\")\n", + " plot_field(TUTORIAL_DIR / \"membrane_20.e\", mu_s, ax, f\"Sample field mu {i}\")\n", + " ax.set_title(f\"Sample field mu {i}\")\n", + "\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As we can see, we set the material parameter to be constant in each element. We can think about this random field of a domain that has some type of defect where the mechanical properties are different.\n", + "\n", + "> Note: This tutorial assumes a working local 4C installation. Please create a symbolic link to your 4C build directory and store it under `/config/4C_build` via this command:\n", + "> ```\n", + "> ln -s /config/4C_build\n", + "> ```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Let's do Monte Carlo integration\n", + "\n", + "from queens.parameters import Parameters\n", + "from queens.data_processors import PvdFile\n", + "from queens.schedulers import Local\n", + "from queens_interfaces.fourc.driver import Fourc\n", + "from queens.models import Simulation\n", + "from queens.iterators import MonteCarlo\n", + "from queens.global_settings import GlobalSettings\n", + "from queens.main import run_iterator\n", + "from queens.utils.io import load_result\n", + "\n", + "# Set the paths for 4C executable and input template file\n", + "fourc_executable = QUEENS_BASE_DIR / \"config\" / \"4C_build\" / \"4C\"\n", + "# input_file_template = TUTORIAL_DIR / \"coarse_plate_dirichlet_template.4C.yaml\"\n", + "input_file_template = QUEENS_BASE_DIR / \"tests\" / \"input_files\" / \"third_party\" / \"fourc\" / \"coarse_plate_dirichlet_template.4C.yaml\"\n", + "input_mesh = TUTORIAL_DIR / \"membrane_20.e\"\n", + "\n", + "if __name__ == \"__main__\":\n", + " with GlobalSettings(\n", + " \"monte_carlo_random_field_4C\", output_dir=\"monte_carlo_random_field_4C\"\n", + " ) as gs:\n", + " parameters = Parameters(MUE=random_field)\n", + "\n", + " # Extract the Cauchy stresses at the last time step\n", + " data_processor = PvdFile(\n", + " field_name=\"element_cauchy_stresses_xyz\",\n", + " file_name_identifier=\"output-structure.pvd\",\n", + " file_options_dict={},\n", + " point_data=False,\n", + " time_steps=[-1],\n", + " )\n", + "\n", + " # How to run 4C\n", + " driver = Fourc(\n", + " parameters=parameters,\n", + " input_templates={\n", + " \"input_file\": input_file_template, # Input file for 4C\n", + " \"material_file\": material_file_template, # File containing the random field\n", + " },\n", + " executable=fourc_executable,\n", + " data_processor=data_processor,\n", + " files_to_copy=[input_mesh], # copy mesh file to experiment directory\n", + " )\n", + "\n", + " # Schedule parallel simulations\n", + " scheduler = Local(num_procs=1, num_jobs=4, experiment_name=gs.experiment_name)\n", + "\n", + " # our model\n", + " model = Simulation(scheduler=scheduler, driver=driver)\n", + " iterator = MonteCarlo(\n", + " seed=1,\n", + " num_samples=100,\n", + " result_description={\"write_results\": True, \"plot_results\": False},\n", + " model=model,\n", + " parameters=parameters,\n", + " global_settings=gs,\n", + " )\n", + "\n", + " run_iterator(iterator, gs)\n", + "\n", + " results = load_result(gs.result_file(\"pickle\"))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Let's plot the mean Cauchy stresses\n", + "\n", + "input_samples = results[\"input_data\"]\n", + "cauchy_stresses = results[\"raw_output_data\"][\"result\"]\n", + "\n", + "\n", + "fig, axes = plt.subplots(1, 3, figsize=(15, 5), subplot_kw={\"projection\": \"3d\"})\n", + "\n", + "plot_field(\n", + " file_path=TUTORIAL_DIR / \"membrane_20.e\",\n", + " field=results[\"mean\"][:, 0],\n", + " ax=axes[0],\n", + " color_bar_title=\"Mean Cauchy Stress xx\"\n", + ")\n", + "\n", + "plot_field(\n", + " file_path=TUTORIAL_DIR / \"membrane_20.e\",\n", + " field=results[\"mean\"][:, 1],\n", + " ax=axes[1],\n", + " color_bar_title=\"Mean Cauchy Stress yy\"\n", + ")\n", + "\n", + "plot_field(\n", + " file_path=TUTORIAL_DIR / \"membrane_20.e\",\n", + " field=results[\"mean\"][:, 3],\n", + " ax=axes[2],\n", + " color_bar_title=\"Mean Cauchy Stress xy\"\n", + ")\n", + "\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Let's have a look at some outputs\n", + "\n", + "fig, axes = plt.subplots(4, 3, figsize=(15, 20), subplot_kw={\"projection\": \"3d\"})\n", + "\n", + "job_ids = [0, 50, 99]\n", + "\n", + "for i, job_id in enumerate(job_ids):\n", + " plot_field(\n", + " file_path=TUTORIAL_DIR / \"membrane_20.e\",\n", + " field=random_field.expanded_representation(input_samples[job_id, :]),\n", + " ax=axes[0, i],\n", + " color_bar_title=f\"Sample {job_id}: Shear modulus\"\n", + " )\n", + "\n", + " plot_field(\n", + " file_path=TUTORIAL_DIR / \"membrane_20.e\",\n", + " field=cauchy_stresses[job_id, :][:, 0],\n", + " ax=axes[1, i],\n", + " color_bar_title=f\"Sample {job_id} Cauchy Stress xx\"\n", + " )\n", + "\n", + " plot_field(\n", + " file_path=TUTORIAL_DIR / \"membrane_20.e\",\n", + " field=cauchy_stresses[job_id, :][:, 1],\n", + " ax=axes[2, i],\n", + " color_bar_title=f\"Sample {job_id} Cauchy Stress yy\"\n", + " )\n", + "\n", + " plot_field(\n", + " file_path=TUTORIAL_DIR / \"membrane_20.e\",\n", + " field=cauchy_stresses[job_id, :][:, 3],\n", + " ax=axes[3, i],\n", + " color_bar_title=f\"Sample {job_id} Cauchy Stress xy\"\n", + " )\n", + "\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We started 4C using QUEENS and prescribed a heterogeneous constitutive parameters! Note that compared to the previous examples, we only had to adapt the parameters! From driver to scheduler to model to iterator, nothing changed, highlighting the modularity of QUEENS." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Gaussian random fields\n", + "\n", + "There are various ways of defining a random field, a common one being Gaussian random field (also known as [Gaussian process](https://en.wikipedia.org/wiki/Gaussian_process))! The main aspect of Gaussian random fields $f(x,y)$, is that at every location $x,y$ of the domain, the random field is normally distributed:\n", + "\n", + "$f(x,y) \\sim \\mathcal{N}(f|\\mu_f(x,y),K((x,y), (x',y')))$\n", + "\n", + "Here $\\mu_f(x,y)$ is the mean value function and $K$ the covariance function. The latter one describes the properties of the random field, such as smoothness, lengthscales, etc.\n", + "\n", + "Gaussian random fields can be sampled using the [Kosambi–Karhunen–Loève](https://en.wikipedia.org/wiki/Kosambi%E2%80%93Karhunen%E2%80%93Lo%C3%A8ve_theorem) theorem:\n", + "\n", + "$$f^{(s)}(x,y) = \\mu_f(x,y)+\\sum_{k=1}^{\\infty} \\sqrt{\\lambda_k} \\xi_k^{(s)} e_k(x,y)$$\n", + "\n", + "where $\\lambda_k$ and $e_k$ are the $k$-th eigenvalues and eigenfunctions of $K$, where $\\boldsymbol{\\xi}^{(s)}$ are samples of Gaussian random variables with zero mean and identity covariance matrix. This approach is similar to [proper orthogonal decomposition](https://en.wikipedia.org/wiki/Proper_orthogonal_decomposition).\n", + "\n", + "For these, we truncate this series at some value $k_{t}$. We select the covariance function\n", + "\n", + "$$K((x,y),(x',y')) = \\sigma_G^2 \\exp\\left(-\\frac{(x-x')^2 + (y-y')^2}{2l^2}\\right)$$\n", + "\n", + "where we set the length scale $l=8.0$ and the output variance $\\sigma_G^2 = 0.03$. The mean function is assumed constant as $\\mu_f(x,y)=0.25$. The number of eigenfunctions $k_t$ is chosen such that the explained variance equals 0.95, i.e. $\\frac{\\sum_{k=1}^{k_t} \\lambda_k}{\\sum_{j=1}^{\\infty} \\lambda_j} \\approx 0.95$.\n", + "\n", + "\n", + "*Enough theory, let's look at some samples:*" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from queens.parameters.random_fields import KarhunenLoeve\n", + "\n", + "np.random.seed(42)\n", + "\n", + "random_field = KarhunenLoeve(\n", + " corr_length=8.0,\n", + " std=0.03,\n", + " mean=0.25,\n", + " explained_variance=0.95,\n", + " coords=mu_rf_parameters,\n", + ")\n", + "\n", + "latent_samples = random_field.draw(3)\n", + "mu_samples = random_field.expanded_representation(latent_samples)\n", + "\n", + "fig, axes = plt.subplots(1, 3, figsize=(15, 5), subplot_kw={\"projection\": \"3d\"})\n", + "\n", + "for i, mu_s in enumerate(mu_samples):\n", + " plot_field(\n", + " file_path=TUTORIAL_DIR / \"membrane_20.e\",\n", + " field=mu_s,\n", + " ax=axes[i],\n", + " color_bar_title=f\"Sample field mu {i}\"\n", + " )\n", + "\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You can see how the samples are distinct from the previous examples. They are more 'wiggly', yet still one can see a correlation structure. This is precisely what the length scale and variance do!\n", + "\n", + "Try it out yourself, decrease/increase the parameters and look at the samples.\n", + "\n", + "What do these parameters do?\n", + "\n", + "
\n", + "\n", + "Answer\n", + "\n", + "According to [The Kernel Cookbook by David Duvenaud](https://www.cs.toronto.edu/~duvenaud/cookbook/):\n", + "\n", + "- The length scale determines the length of the 'wiggles' in your function. In general, you won't be able to extrapolate more than $l$ units away from your data.\n", + "- The output variance $\\sigma_G^2$ determines the average distance of your function away from its mean. Every kernel has this parameter out in front; it's just a scale factor.\n", + "
\n", + "\n", + "Look at the cutoff value $k_t$ (hint: `random_field.dimension`). How does this change depending on the length scale?\n", + "\n", + "
\n", + "\n", + "Answer\n", + "\n", + "With smaller length scales, the 'wiggle'-frequency increases, hence more components are needed to explain the desired variance in eigenvalues of 95%.\n", + "
\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Nice, now let us repeat the Monte Carlo experiment using the Gaussian random field described above!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "if __name__ == \"__main__\":\n", + " with GlobalSettings(\n", + " \"monte_carlo_gaussian_random_field_4C\", output_dir=\"monte_carlo_gaussian_random_field_4C\"\n", + " ) as gs:\n", + " parameters = Parameters(MUE=random_field)\n", + "\n", + " # Extract the Cauchy stresses at the last time step\n", + " data_processor = PvdFile(\n", + " field_name=\"element_cauchy_stresses_xyz\",\n", + " file_name_identifier=\"output-structure.pvd\",\n", + " file_options_dict={},\n", + " point_data=False,\n", + " time_steps=[-1],\n", + " files_to_be_deleted_regex_lst=[\n", + " \"output.control\",\n", + " \"*.mesh.s0\",\n", + " \"output-vtk-files/*\",\n", + " ], # These files are deleted\n", + " )\n", + "\n", + " # How to run 4C\n", + " driver = Fourc(\n", + " parameters=parameters,\n", + " input_templates={\n", + " \"input_file\": input_file_template, # Input file for 4C\n", + " \"material_file\": material_file_template, # File containing the random field\n", + " },\n", + " executable=fourc_executable,\n", + " data_processor=data_processor,\n", + " files_to_copy=[input_mesh], # copy mesh file to experiment directory\n", + " )\n", + "\n", + " # Schedule parallel simulations\n", + " scheduler = Local(num_procs=1, num_jobs=4, experiment_name=gs.experiment_name)\n", + "\n", + " # Our model\n", + " model = Simulation(scheduler=scheduler, driver=driver)\n", + " iterator = MonteCarlo(\n", + " seed=1,\n", + " num_samples=100,\n", + " result_description={\"write_results\": True, \"plot_results\": False},\n", + " model=model,\n", + " parameters=parameters,\n", + " global_settings=gs,\n", + " )\n", + "\n", + " run_iterator(iterator, gs)\n", + "\n", + " results = load_result(gs.result_file(\"pickle\"))\n", + "\n", + "# Let's plot the mean Cauchy stresses\n", + "\n", + "input_samples = results[\"input_data\"]\n", + "\n", + "fig, axes = plt.subplots(1, 3, figsize=(15, 5), subplot_kw={\"projection\": \"3d\"})\n", + "\n", + "plot_field(\n", + " file_path=TUTORIAL_DIR / \"membrane_20.e\",\n", + " field=results[\"mean\"][:, 0],\n", + " ax=axes[0],\n", + " color_bar_title=\"Mean Cauchy Stress xx\"\n", + ")\n", + "\n", + "plot_field(\n", + " file_path=TUTORIAL_DIR / \"membrane_20.e\",\n", + " field=results[\"mean\"][:, 1],\n", + " ax=axes[1],\n", + " color_bar_title=\"Mean Cauchy Stress yy\"\n", + ")\n", + "\n", + "plot_field(\n", + " file_path=TUTORIAL_DIR / \"membrane_20.e\",\n", + " field=results[\"mean\"][:, 3],\n", + " ax=axes[2],\n", + " color_bar_title=\"Mean Cauchy Stress xy\"\n", + ")\n", + "\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Why QUEENS?\n", + "\n", + "- Hide complexity in handling complex simulation codes (same as in the last example)\n", + "- QUEENS hides also the complexity of working with random fields:\n", + " - Monte Carlo iterator does not care that it is a random field\n", + " - We can still use the native 4C driver even though we now have 2 input files (one for the simulation parameters, one for the random field)\n", + " \n", + "## Let's play around\n", + "\n", + "Let your creativity flow and try out stuff.\n", + "\n", + "### Inspirations\n", + "- Create your own random field definition, maybe a sine wave? (Keep in mind the parameter has to remain positive!)\n", + "- Take a look at some simulation output of the Monte Carlo run (change the data processor attribute `files_to_be_deleted_regex_lst` to `None`)\n", + "- Rerun the Monte Carlo analysis with different length scales\n", + "- Rerun the Monte Carlo analysis with different numbers of samples\n", + "- Plot the variance estimate from the Monte Carlo run\n", + "- Change the length scale and look at the Monte Carlo mean value for multiple runs with the same number of samples" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "psdouu", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.0" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/tutorials/4_quantifying_uncertainty_due_to_heterogeneous_material_fields/coarse_plate_dirichlet_template.4C.yaml b/tutorials/4_quantifying_uncertainty_due_to_heterogeneous_material_fields/coarse_plate_dirichlet_template.4C.yaml new file mode 100644 index 000000000..f42be5078 --- /dev/null +++ b/tutorials/4_quantifying_uncertainty_due_to_heterogeneous_material_fields/coarse_plate_dirichlet_template.4C.yaml @@ -0,0 +1,100 @@ +TITLE: 4C Random field job {{ job_id }} +PROBLEM TYPE: + PROBLEMTYPE: "Structure" +IO: + OUTPUT_SPRING: true + STRUCT_STRESS: Cauchy +IO/RUNTIME VTK OUTPUT: + INTERVAL_STEPS: 1 +IO/RUNTIME VTK OUTPUT/STRUCTURE: + OUTPUT_STRUCTURE: true + DISPLACEMENT: true + STRESS_STRAIN: true +STRUCTURAL DYNAMIC: + INT_STRATEGY: "Standard" + RESEVERYERGY: 1 + RESTARTEVERY: 0 + TIMESTEP: 0.025 + NUMSTEP: 40 + MAXTIME: 1 + M_DAMP: 0.001 + K_DAMP: 0.001 + TOLDISP: 1e-06 + TOLRES: 1e-06 + MAXITER: 200 + PREDICT: "ConstDisVelAcc" + LINEAR_SOLVER: 1 +SOLVER 1: + SOLVER: "UMFPACK" + NAME: "direct solver" +MATERIALS: + - MAT: 10 + MAT_Membrane_ElastHyper: + NUMMAT: 1 + MATIDS: [11] + DENS: 1e-09 + - MAT: 11 + ELAST_IsoNeoHooke: + MUE: + from_file: "{{ material_file }}" +FUNCT1: + - COMPONENT: 0 + SYMBOLIC_FUNCTION_OF_SPACE_TIME: "l" + - VARIABLE: 0 + NAME: "l" + TYPE: "linearinterpolation" + NUMPOINTS: 2 + TIMES: [0, 1] + VALUES: [0, 2.5] +FUNCT2: + - COMPONENT: 0 + SYMBOLIC_FUNCTION_OF_SPACE_TIME: "l" + - VARIABLE: 0 + NAME: "l" + TYPE: "linearinterpolation" + NUMPOINTS: 2 + TIMES: [0, 1] + VALUES: [0, 2.5] +DESIGN LINE DIRICH CONDITIONS: + - E: 3 + ENTITY_TYPE: node_set_id + NUMDOF: 3 + ONOFF: [0, 1, 1] + VAL: [0, 0, 0] + FUNCT: [0, 0, 0] + - E: 2 + ENTITY_TYPE: node_set_id + NUMDOF: 3 + ONOFF: [1, 0, 1] + VAL: [0, 0, 0] + FUNCT: [0, 0, 0] + - E: 1 + ENTITY_TYPE: node_set_id + NUMDOF: 3 + ONOFF: [0, 1, 1] + VAL: [0, 1, 0] + FUNCT: [0, 2, 0] + - E: 4 + ENTITY_TYPE: node_set_id + NUMDOF: 3 + ONOFF: [1, 0, 1] + VAL: [1, 0, 0] + FUNCT: [1, 0, 0] +DESIGN SURF DIRICH CONDITIONS: + - E: 1 + ENTITY_TYPE: element_block_id + NUMDOF: 3 + ONOFF: [0, 0, 1] + VAL: [0, 0, 0] + FUNCT: [0, 0, 0] +STRUCTURE GEOMETRY: + FILE: ../membrane_20.e + SHOW_INFO: none + ELEMENT_BLOCKS: + - ID: 1 + MEMBRANE4: + QUAD4: + MAT: 10 + KINEM: nonlinear + THICK: 2.5 + STRESS_STRAIN: plane_stress diff --git a/tutorials/4_quantifying_uncertainty_due_to_heterogeneous_material_fields/membrane_20.e b/tutorials/4_quantifying_uncertainty_due_to_heterogeneous_material_fields/membrane_20.e new file mode 100644 index 000000000..6deaac7a1 Binary files /dev/null and b/tutorials/4_quantifying_uncertainty_due_to_heterogeneous_material_fields/membrane_20.e differ diff --git a/tutorials/4_quantifying_uncertainty_due_to_heterogeneous_material_fields/plot_input_random_field.py b/tutorials/4_quantifying_uncertainty_due_to_heterogeneous_material_fields/plot_input_random_field.py new file mode 100644 index 000000000..0c7cfa48b --- /dev/null +++ b/tutorials/4_quantifying_uncertainty_due_to_heterogeneous_material_fields/plot_input_random_field.py @@ -0,0 +1,55 @@ +# +# SPDX-License-Identifier: LGPL-3.0-or-later +# Copyright (c) 2024-2025, QUEENS contributors. +# +# This file is part of QUEENS. +# +# QUEENS is free software: you can redistribute it and/or modify it under the terms of the GNU +# Lesser General Public License as published by the Free Software Foundation, either version 3 of +# the License, or (at your option) any later version. QUEENS is distributed in the hope that it will +# be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You +# should have received a copy of the GNU Lesser General Public License along with QUEENS. If not, +# see . +# +"""Plotting functions for tutorial 4.""" +import matplotlib.pyplot as plt +import mpl_toolkits +import numpy as np +import pyvista as pv +from mpl_toolkits.mplot3d.art3d import Poly3DCollection +from mpl_toolkits.mplot3d.axes3d import Axes3D + + +def plot_field(file_path: str, field: np.ndarray, ax: Axes3D, color_bar_title: str) -> None: + """Plot cell-wise field values on the surface of a mesh. + + Args: + file_path: Path to the mesh file readable by PyVista. + field: Cell data array of shape ``(num_cells,)``. + ax: 3D Matplotlib axis to plot into. + color_bar_title: Label shown next to the color bar. + """ + mesh = pv.read(file_path)[0][0] + mesh.cell_data["field"] = field + cell_values = field + + faces = mesh.extract_surface().faces.reshape((-1, 5))[:, 1:] + vertices = mesh.points + + cmap = plt.get_cmap("viridis") + norm = plt.Normalize(vmin=cell_values.min(), vmax=cell_values.max()) + colors = [cmap(norm(val)) for val in cell_values] + + poly3d = [[vertices[idx] for idx in face] for face in faces] + collection = Poly3DCollection(poly3d, facecolors=colors, edgecolor="k", linewidths=0.2) + ax.add_collection3d(collection) + + ax.auto_scale_xyz(vertices[:, 0], vertices[:, 1], vertices[:, 2]) + ax.set_axis_off() + + sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm) + sm.set_array([]) + plt.colorbar(sm, ax=ax, fraction=0.02, pad=0.04, label=color_bar_title) + + ax.view_init(elev=90, azim=-90) diff --git a/tutorials/4_quantifying_uncertainty_due_to_heterogeneous_material_fields/random_field.py b/tutorials/4_quantifying_uncertainty_due_to_heterogeneous_material_fields/random_field.py new file mode 100644 index 000000000..3106223a5 --- /dev/null +++ b/tutorials/4_quantifying_uncertainty_due_to_heterogeneous_material_fields/random_field.py @@ -0,0 +1,90 @@ +# +# SPDX-License-Identifier: LGPL-3.0-or-later +# Copyright (c) 2024-2025, QUEENS contributors. +# +# This file is part of QUEENS. +# +# QUEENS is free software: you can redistribute it and/or modify it under the terms of the GNU +# Lesser General Public License as published by the Free Software Foundation, either version 3 of +# the License, or (at your option) any later version. QUEENS is distributed in the hope that it will +# be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You +# should have received a copy of the GNU Lesser General Public License along with QUEENS. If not, +# see . +# +"""CustomRandomField module for tutorial 4.""" +from collections.abc import Callable + +import numpy as np + +from queens.distributions._distribution import Continuous, Discrete +from queens.parameters.random_fields._random_field import RandomField + + +class CustomRandomField(RandomField): + """Random field defined by an explicit expansion function. + + Attributes: + latent_distribution: QUEENS distribution object of random field + expansion: Callable mapping ``(theta, coords) -> field``. + coordinates: Coordinates at which the random field is evaluated. + """ + + def __init__( + self, coords: dict, latent_distribution: Continuous | Discrete, expansion: Callable + ) -> None: + """Initialize the custom random field. + + Args: + coords: Dictionary with coordinates at which the random field is evaluated. + latent_distribution: Distribution of the latent parameters ``theta``. + expansion: Function computing field values from ``(theta, coords)``. + """ + super().__init__(coords, latent_distribution, latent_distribution.dimension) + self.latent_distribution = latent_distribution + self.expansion = expansion + self.coordinates = coords["coords"] + + def draw(self, num_samples: int) -> np.ndarray: + """Draw ``num_samples`` samples of the latent space. + + Args: + num_samples (int): Batch size of samples to draw + """ + return self.latent_distribution.draw(num_samples) + + def expanded_representation(self, samples: np.ndarray) -> np.ndarray: + """Expand the random field realization. + + Args: + samples (np.array): Latent space variables to be expanded into a random field + """ + if samples.ndim == 1: + return self.expansion(samples, self.coordinates) + + expansions = [] + + for sample in samples: + expansions.append(self.expansion(sample, self.coordinates)) + + return np.array(expansions) + + def logpdf(self, samples: np.ndarray) -> np.ndarray: + """Get joint logpdf of latent space. + + Args: + samples (np.array): Sample to evaluate logpdf + """ + return self.latent_distribution.logpdf(samples) + + def grad_logpdf(self, samples: np.ndarray) -> np.ndarray: + """Get gradient of joint logpdf of latent space. + + Args: + samples (np.array): Sample to evaluate gradient of logpdf + """ + raise NotImplementedError() + + def latent_gradient(self, upstream_gradient: np.ndarray) -> np.ndarray: + """Gradient of the field with respect to the latent variables.""" + raise NotImplementedError() diff --git a/tutorials/__init__.py b/tutorials/__init__.py new file mode 100644 index 000000000..32b0ac211 --- /dev/null +++ b/tutorials/__init__.py @@ -0,0 +1,14 @@ +# +# SPDX-License-Identifier: LGPL-3.0-or-later +# Copyright (c) 2024-2025, QUEENS contributors. +# +# This file is part of QUEENS. +# +# QUEENS is free software: you can redistribute it and/or modify it under the terms of the GNU +# Lesser General Public License as published by the Free Software Foundation, either version 3 of +# the License, or (at your option) any later version. QUEENS is distributed in the hope that it will +# be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You +# should have received a copy of the GNU Lesser General Public License along with QUEENS. If not, +# see . +# diff --git a/tutorials/utils.py b/tutorials/utils.py new file mode 100644 index 000000000..da9f4cc10 --- /dev/null +++ b/tutorials/utils.py @@ -0,0 +1,31 @@ +# +# SPDX-License-Identifier: LGPL-3.0-or-later +# Copyright (c) 2024-2025, QUEENS contributors. +# +# This file is part of QUEENS. +# +# QUEENS is free software: you can redistribute it and/or modify it under the terms of the GNU +# Lesser General Public License as published by the Free Software Foundation, either version 3 of +# the License, or (at your option) any later version. QUEENS is distributed in the hope that it will +# be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You +# should have received a copy of the GNU Lesser General Public License along with QUEENS. If not, +# see . +# +"""Utility functions for tutorials.""" +from pathlib import Path + + +def find_repo_root(start: Path) -> Path: + """Find the repository root by going upwards from a starting directory. + + Args: + start (Path): Starting directory to begin the upward search. + + Returns: + Path: Path to the detected repository root (containing ``pyproject.toml``). + """ + for p in [start, *start.parents]: + if (p / "pyproject.toml").exists(): + return p + raise RuntimeError(f"Could not find repo root from {start}")