Add new tutorial: two-scale-heat-conduction (#343)

* Initial commit - add Nutils macro and micro participants * Initial code for macro and micro cases with DuMux * Add cmake file for compiling DuMux code * Remove DuMux case for now * Initial content in README and simplifying config scripts * Use standard name of preCICE config * Add cleaning script * Add clean scripts for macro- and micro- Nutils * Use copies of Nutils objects while getting the state * Add image and modify clean-tutorial.sh * Reducing computational time of the macro-micro case and adding content to README * Making the case cheaper and finalizing README * Formatting * Update two-scale-heat-conduction/README.md Co-authored-by: Gerasimos Chourdakis <[email protected]> * Apply suggestions from code review Co-authored-by: Gerasimos Chourdakis <[email protected]> * Add run scripts for macro and micro Nutils participants * Simplifying README and formatting precice-config * Add tutorial prefix name to image filename * Rename image file * Add shebang in clean-tutorial.sh * Use quotes while passing numbers of procs to run.sh * Incorporating changes * Add correct file relative paths * Port case to v3 * Add results image in README and do some formatting * Add results image * Add micro simulation results and descriptions of outputs * Change width of images in README * Remove width parameter, it does not work * Small typos --------- Co-authored-by: Gerasimos Chourdakis <[email protected]>
precice · Sep 5, 2023 · 847a6f0 · 847a6f0
1 parent f67ae5a
commit 847a6f0
Show file tree

Hide file tree

Showing 14 changed files with 698 additions and 0 deletions.
diff --git a/two-scale-heat-conduction/README.md b/two-scale-heat-conduction/README.md
@@ -0,0 +1,63 @@
+---
+title: Two-scale heat conduction
+permalink: tutorials-two-scale-heat-conduction.html
+keywords: Macro-micro, Micro Manager, Nutils, Heat conduction
+summary: We solve a two-scale heat conduction problem with a predefined micro structure of two materials. One macro simulation is coupled to several micro simulations using the Micro Manager.
+---
+
+{% note %}
+Get the [case files of this tutorial](https://github.com/precice/tutorials/tree/master/two-scale-heat-conduction). Read how in the [tutorials introduction](https://www.precice.org/tutorials.html).
+{% endnote %}
+
+## Setup
+
+This tutorial solves a heat conduction problem on a 2D domain which has an underlying micro-structure. This micro-structure changes the constituent quantities necessary for solving the problem on the macro scale. This leads to a two-scale problem with one macro-scale simulation and several micro-scale simulations.
+
+![Case setup of two-scale-heat-conduction case](images/tutorials-two-scale-heat-conduction-macro-micro-schematic.png)
+
+At each Gauss point of the macro domain there exists a micro simulation. The macro problem is one participant, which is coupled to many micro simulations. Each micro simulation is not an individual coupling participant, instead we use a managing software which controls all the micro simulations and their coupling via preCICE. The case is chosen from the first example case in the publication
+
+*Bastidas, Manuela & Bringedal, Carina & Pop, Iuliu Sorin (2021), A two-scale iterative scheme for a phase-field model for precipitation and dissolution in porous media. Applied Mathematics and Computation. 396. 125933. [10.1016/j.amc.2020.125933](https://doi.org/10.1016/j.amc.2020.125933)*.
+
+## Available solvers and dependencies
+
+* Both the macro and micro simulations are solved using the finite element library [Nutils](https://nutils.org/install.html) v7.
+* The macro simulation is written in Python, so it requires the [Python bindings of preCICE](https://precice.org/installation-bindings-python.html)
+* The [Micro Manager](https://precice.org/tooling-micro-manager-installation.html) controls all micro-simulations and facilitates coupling via preCICE. Use the [develop](https://github.com/precice/micro-manager/tree/develop) branch of the Micro Manager.
+
+## Running the simulation
+
+Run the macro problem:
+
+```bash
+cd macro-nutils
+./run.sh
+```
+
+Check the Micro Manager [configuration](https://precice.org/tooling-micro-manager-configuration.html) and [running](https://precice.org/tooling-micro-manager-running.html) documentation to understand how to set it up and launch it. There is a Python script `run-micro-problems.py` in the tutorial directory to directly run the Micro Manager. This script imports the Micro Manager, and calls its `initialize()` and `solve()` methods. The Micro Manager can be run via this script in serial or parallel. Run it in serial:
+
+```bash
+cd micro-nutils
+./run.sh -s
+```
+
+Run it in parallel:
+
+```bash
+cd micro-nutils
+./run.sh -p <num_procs>
+```
+
+The `num_procs` needs to fit the decomposition specified in the `micro-manager-config.json` (default: two ranks).
+
+Even though the case setup and involved physics is simple, each micro simulation is an instance of Nutils, which usually has a moderately high computation time. If the Micro Manager is run on 2 processors, the total runtime is approximately 25-30 minutes. Do not run the Micro Manager in serial, because the runtime will be several hours.
+
+## Post-processing
+
+![Results of two-scale-heat-conduction case](images/tutorials-two-scale-heat-conduction-results.png)
+
+The participant `macro-nutils` outputs `macro-*.vtk` files which can be viewed in ParaView to see the macro concentration field. The Micro Manager uses the [export functionality](https://precice.org/configuration-export.html#enabling-exporters) of preCICE to output micro simulation data and [adaptivity related data](https://precice.org/tooling-micro-manager-configuration.html#adding-adaptivity-in-the-precice-xml-configuration) to VTU files which can be viewed in ParaView. To view the data on each micro simulation, create a Glyph on the Micro Manager VTU data. In the figure above, micro-scale porosity is shown. For a lower concentration value, the porosity increases (in the lower left corner).
+
+![Evolving micro simulations](images/tutorials-two-scale-heat-conduction-evolving-micro-simulations.png)
+
+The micro simulations themselves have a circular micro structure which is resolved in every time step. To output VTK files for each micro simulation, uncomment the `output()` function in the file `micro-nutils/micro.py`. The figure above shows the changing phase field used to represent the circular micro structure and the diffuse interface width.
diff --git a/two-scale-heat-conduction/clean-tutorial.sh b/two-scale-heat-conduction/clean-tutorial.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+../tools/clean-tutorial-base.sh
diff --git a/...ction/images/tutorials-two-scale-heat-conduction-evolving-micro-simulations.png b/...ction/images/tutorials-two-scale-heat-conduction-evolving-micro-simulations.png
diff --git a/...conduction/images/tutorials-two-scale-heat-conduction-macro-micro-schematic.png b/...conduction/images/tutorials-two-scale-heat-conduction-macro-micro-schematic.png
diff --git a/two-scale-heat-conduction/images/tutorials-two-scale-heat-conduction-results.png b/two-scale-heat-conduction/images/tutorials-two-scale-heat-conduction-results.png
diff --git a/two-scale-heat-conduction/macro-nutils/clean.sh b/two-scale-heat-conduction/macro-nutils/clean.sh
@@ -0,0 +1,6 @@
+#!/bin/sh
+set -e -u
+
+. ../../tools/cleaning-tools.sh
+
+clean_nutils .
diff --git a/two-scale-heat-conduction/macro-nutils/macro.py b/two-scale-heat-conduction/macro-nutils/macro.py
@@ -0,0 +1,166 @@
+#! /usr/bin/env python3
+#
+# Heat conduction problem with two materials.
+
+from nutils import mesh, function, solver, export, cli
+import treelog
+import numpy as np
+import precice
+
+
+def main():
+    """
+    2D heat conduction problem solved on a rectangular domain.
+    The material consists of a mixture of two materials "g" and "s".
+    """
+    is_coupled_case = True  # If False, single-physics problem is solved
+
+    # Original case from Bastidas et al.
+    # topo, geom = mesh.rectilinear([np.linspace(0, 1.0, 9), np.linspace(0, 0.5, 5)])
+    topo, geom = mesh.rectilinear([np.linspace(0, 1.0, 5), np.linspace(0, 0.5, 4)])
+
+    ns = function.Namespace(fallback_length=2)
+    ns.x = geom
+    ns.basis = topo.basis('std', degree=1)
+    ns.u = 'basis_n ?solu_n'  # the field to be solved for
+
+    ns.phi = 'basis_n ?solphi_n'  # porosity, to be read from preCICE
+    phi = 0.5  # initial value of porosity
+
+    ns.kbasis = topo.basis('std', degree=1).vector(topo.ndims).vector(topo.ndims)
+    # conductivity matrix, to be read from preCICE (read as two vectors and then patched to form a matrix)
+    ns.k_ij = 'kbasis_nij ?solk_n'
+    k = 1.0  # initial value of conductivity
+
+    ns.rhos = 1.0  # density of material s
+    ns.rhog = 1.0  # density of material g
+    ns.dudt = 'basis_n (?solu_n - ?solu0_n) / ?dt'  # implicit Euler time stepping formulation
+
+    # Dirichlet boundary condition at lower left corner
+    ns.usource = 0.0
+
+    # Initial condition in the domain
+    ns.uinitial = 0.5
+
+    if is_coupled_case:
+        participant = precice.Participant("macro-heat", "../precice-config.xml", 0, 1)
+        mesh_name = "macro-mesh"
+
+        # Define Gauss points on entire domain as coupling mesh (volume coupling from macro side)
+        couplingsample = topo.sample('gauss', degree=2)  # mesh vertices are Gauss points
+        vertex_ids = participant.set_mesh_vertices(mesh_name, couplingsample.eval(ns.x))
+    else:
+        sqrphi = topo.integral((ns.phi - phi) ** 2, degree=1)
+        solphi = solver.optimize('solphi', sqrphi, droptol=1E-12)
+
+        sqrk = topo.integral(((ns.k - k * np.eye(2)) * (ns.k - k * np.eye(2))).sum([0, 1]), degree=2)
+        solk = solver.optimize('solk', sqrk, droptol=1E-12)
+
+    # Define the weak form of the heat conduction problem with two materials
+    res = topo.integral('((rhos phi + (1 - phi) rhog) basis_n dudt + k_ij basis_n,i u_,j) d:x' @ ns, degree=2)
+
+    # Set Dirichlet boundary conditions
+    sqr = topo.boundary['left'].boundary['bottom'].integral('(u - usource)^2 d:x' @ ns, degree=2)
+    cons = solver.optimize('solu', sqr, droptol=1e-15)
+
+    # Set domain to initial condition
+    sqr = topo.integral('(u - uinitial)^2' @ ns, degree=2)
+    solu0 = solver.optimize('solu', sqr)
+
+    if is_coupled_case:
+        if participant.requires_initial_data():
+            concentrations = couplingsample.eval('u' @ ns, solu=solu0)
+            participant.write_data(mesh_name, "concentration", vertex_ids, concentrations)
+
+        participant.initialize()
+        dt = participant.get_max_time_step_size()
+    else:
+        dt = 1.0E-2
+
+    ns.dt = dt
+    n = n_checkpoint = 0
+    t = t_checkpoint = 0
+    t_out = 0.01
+    t_end = 0.25  # Only relevant when single physics case is run
+    n_out = int(t_out / dt)
+    n_t = int(t_end / dt)
+
+    # Prepare the post processing sample
+    bezier = topo.sample('bezier', 2)
+
+    # Output of initial state
+    x, u = bezier.eval(['x_i', 'u'] @ ns, solu=solu0)
+    with treelog.add(treelog.DataLog()):
+        export.vtk('macro-0', bezier.tri, x, u=u)
+
+    if is_coupled_case:
+        is_coupling_ongoing = participant.is_coupling_ongoing()
+    else:
+        is_coupling_ongoing = True
+
+    while is_coupling_ongoing:
+        if is_coupled_case:
+            if participant.requires_writing_checkpoint():
+                solu_checkpoint = solu0
+                t_checkpoint = t
+                n_checkpoint = n
+
+            # Read porosity and apply it to the existing solution
+            poro_data = participant.read_data(mesh_name, "porosity", vertex_ids, dt)
+            poro_coupledata = couplingsample.asfunction(poro_data)
+            sqrphi = couplingsample.integral((ns.phi - poro_coupledata) ** 2)
+            solphi = solver.optimize('solphi', sqrphi, droptol=1E-12)
+
+            # Read conductivity and apply it to the existing solution
+            k_00_c = couplingsample.asfunction(participant.read_data(mesh_name, "k_00", vertex_ids, dt))
+            k_01_c = couplingsample.asfunction(participant.read_data(mesh_name, "k_01", vertex_ids, dt))
+            k_10_c = couplingsample.asfunction(participant.read_data(mesh_name, "k_10", vertex_ids, dt))
+            k_11_c = couplingsample.asfunction(participant.read_data(mesh_name, "k_11", vertex_ids, dt))
+
+            conductivity = function.asarray([[k_00_c, k_01_c], [k_10_c, k_11_c]])
+            sqrk = couplingsample.integral(((ns.k - conductivity) * (ns.k - conductivity)).sum([0, 1]))
+            solk = solver.optimize('solk', sqrk, droptol=1E-12)
+
+        solu = solver.solve_linear('solu', res, constrain=cons,
+                                   arguments=dict(solu0=solu0, dt=dt, solphi=solphi, solk=solk))
+
+        if is_coupled_case:
+            # Collect values of field u and write them to preCICE
+            concentration = couplingsample.eval('u' @ ns, solu=solu)
+            participant.write_data(mesh_name, "concentration", vertex_ids, concentration)
+
+            participant.advance(dt)
+            precice_dt = participant.get_max_time_step_size()
+            dt = min(precice_dt, dt)
+
+        n += 1
+        t += dt
+        solu0 = solu
+
+        if is_coupled_case:
+            is_coupling_ongoing = participant.is_coupling_ongoing()
+
+            if participant.requires_reading_checkpoint():
+                solu0 = solu_checkpoint
+                t = t_checkpoint
+                n = n_checkpoint
+            else:
+                if n % n_out == 0:
+                    x, phi, k, u = bezier.eval(['x_i', 'phi', 'k', 'u'] @ ns, solphi=solphi, solk=solk, solu=solu)
+                    with treelog.add(treelog.DataLog()):
+                        export.vtk('macro-' + str(n), bezier.tri, x, u=u, phi=phi, K=k)
+        else:
+            if n % n_out == 0:
+                x, phi, k, u = bezier.eval(['x_i', 'phi', 'k', 'u'] @ ns, solphi=solphi, solk=solk, solu=solu)
+                with treelog.add(treelog.DataLog()):
+                    export.vtk('macro-' + str(n), bezier.tri, x, u=u, phi=phi, K=k)
+
+            if n >= n_t:
+                is_coupling_ongoing = False
+
+    if is_coupled_case:
+        participant.finalize()
+
+
+if __name__ == '__main__':
+    cli.run(main)
diff --git a/two-scale-heat-conduction/macro-nutils/run.sh b/two-scale-heat-conduction/macro-nutils/run.sh
@@ -0,0 +1,4 @@
+#!/bin/sh
+set -e -u
+
+NUTILS_RICHOUTPUT=no python3 macro.py
diff --git a/two-scale-heat-conduction/micro-nutils/clean.sh b/two-scale-heat-conduction/micro-nutils/clean.sh
@@ -0,0 +1,6 @@
+#!/bin/sh
+set -e -u
+
+. ../../tools/cleaning-tools.sh
+
+clean_nutils .
diff --git a/two-scale-heat-conduction/micro-nutils/micro-manager-config.json b/two-scale-heat-conduction/micro-nutils/micro-manager-config.json
@@ -0,0 +1,26 @@
+{
+    "micro_file_name": "micro",
+    "coupling_params": {
+            "participant_name": "Micro-Manager",
+            "config_file_name": "../precice-config.xml",
+            "macro_mesh_name": "macro-mesh",
+            "write_data_names": {"k_00": "scalar", "k_01": "scalar", "k_10": "scalar", "k_11": "scalar", "porosity": "scalar"},
+            "read_data_names": {"concentration": "scalar"}
+    },
+    "simulation_params": {
+      "macro_domain_bounds": [0.0, 1.0, 0.0, 0.5],
+      "decomposition": [2, 1],
+      "adaptivity": {
+        "type": "global",
+        "data": ["k_00", "k_11", "porosity", "concentration"],
+        "history_param": 0.1,
+        "coarsening_constant": 0.2,
+        "refining_constant": 0.05,
+        "every_implicit_iteration": "False",
+        "similarity_measure": "L2rel"
+      }
+    },
+    "diagnostics": {
+      "data_from_micro_sims": {"grain_size": "scalar"}
+    }
+}