diff --git a/two-scale-heat-conduction/README.md b/two-scale-heat-conduction/README.md new file mode 100644 index 000000000..5f6ef033d --- /dev/null +++ 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 +``` + +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 new file mode 100755 index 000000000..f60ad945b --- /dev/null +++ b/two-scale-heat-conduction/clean-tutorial.sh @@ -0,0 +1,2 @@ +#!/bin/bash +../tools/clean-tutorial-base.sh diff --git a/two-scale-heat-conduction/images/tutorials-two-scale-heat-conduction-evolving-micro-simulations.png b/two-scale-heat-conduction/images/tutorials-two-scale-heat-conduction-evolving-micro-simulations.png new file mode 100644 index 000000000..602b4f6f5 Binary files /dev/null and b/two-scale-heat-conduction/images/tutorials-two-scale-heat-conduction-evolving-micro-simulations.png differ diff --git a/two-scale-heat-conduction/images/tutorials-two-scale-heat-conduction-macro-micro-schematic.png b/two-scale-heat-conduction/images/tutorials-two-scale-heat-conduction-macro-micro-schematic.png new file mode 100644 index 000000000..f6a8a538e Binary files /dev/null and b/two-scale-heat-conduction/images/tutorials-two-scale-heat-conduction-macro-micro-schematic.png differ 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 new file mode 100644 index 000000000..ce00f5c6a Binary files /dev/null and b/two-scale-heat-conduction/images/tutorials-two-scale-heat-conduction-results.png differ diff --git a/two-scale-heat-conduction/macro-nutils/clean.sh b/two-scale-heat-conduction/macro-nutils/clean.sh new file mode 100755 index 000000000..6893c1ea5 --- /dev/null +++ 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 new file mode 100644 index 000000000..56f34197e --- /dev/null +++ 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 new file mode 100755 index 000000000..b646e5ab0 --- /dev/null +++ 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 new file mode 100755 index 000000000..6893c1ea5 --- /dev/null +++ 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 new file mode 100644 index 000000000..74c1a4d4f --- /dev/null +++ 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"} + } +} diff --git a/two-scale-heat-conduction/micro-nutils/micro.py b/two-scale-heat-conduction/micro-nutils/micro.py new file mode 100644 index 000000000..7097a5781 --- /dev/null +++ b/two-scale-heat-conduction/micro-nutils/micro.py @@ -0,0 +1,286 @@ +""" +Micro simulation +In this script we solve the Laplace equation with a grain depicted by a phase field on a square domain :math:`Ω` +with boundary :math:`Γ`, subject to periodic boundary conditions in both dimensions +""" +import math + +from nutils import mesh, function, solver, export, cli +import treelog +import numpy as np +from copy import deepcopy + + +class MicroSimulation: + + def __init__(self): + """ + Constructor of MicroSimulation class. + """ + # Initial parameters + + # self._nelems = 10 # Elements in one direction (original case from Bastidas et al.) + self._nelems = 6 # Elements in one direction + + self._ref_level = 3 # Number of levels of mesh refinement + self._r_initial = 0.4 # Initial radius of the grain + + # Interpolation order for phi and u + self._degree_phi = 2 + self._degree_u = 2 + + # Set up mesh with periodicity in both X and Y directions + self._topo, self._geom = mesh.rectilinear([np.linspace(-0.5, 0.5, self._nelems + 1)] * 2, periodic=(0, 1)) + self._topo_coarse = self._topo # Save original coarse topology to use to re-refinement + + self._solu = None # Solution of weights for which cell problem is solved for + self._solphi = None # Solution of phase field + self._solphinm1 = None # Solution of phase field at t_{n-1} + self._solphi_checkpoint = None # Save the current solution of the phase field as a checkpoint. + self._topo_checkpoint = None # Save the refined mesh as a checkpoint. + self._ucons = None + self._first_iter_done = False + self._initial_condition_is_set = False + self._k_nm1 = None # Average effective conductivity of last time step + + # Define initial namespace + self._ns = function.Namespace() + self._ns.x = self._geom + + self._ns.phibasis = self._topo.basis('std', degree=self._degree_phi) + self._ns.phi = 'phibasis_n ?solphi_n' # Initial phase field + self._ns.coarsephibasis = self._topo_coarse.basis('std', degree=self._degree_phi) + self._ns.coarsephi = 'coarsephibasis_n ?coarsesolphi_n' # Phase field on original coarse topology + self._ns.lam = (4 / self._nelems) / (2 ** self._ref_level) + self._ns.coarselam = 3 / self._nelems + + # Initialize phase field + solphi = self._get_analytical_phasefield( + self._topo, + self._ns, + self._degree_phi, + self._ns.coarselam, + self._r_initial) + + # Refine the mesh + self._topo, self._solphi = self._refine_mesh(self._topo, solphi) + self._reinitialize_namespace(self._topo) + self._initial_condition_is_set = True + + # Initialize phase field once more on refined topology + solphi = self._get_analytical_phasefield(self._topo, self._ns, self._degree_phi, self._ns.lam, self._r_initial) + + target_porosity = 1 - math.pi * self._r_initial ** 2 + print("Target amount of void space = {}".format(target_porosity)) + + dt_initial = 1E-3 + + # Solve Allen-Cahn equation till we reach target porosity value + psi = 0 + while psi < target_porosity: + print("Solving Allen-Cahn equation to achieve initial target grain structure") + solphi = self._solve_allen_cahn(self._topo, solphi, 0.5, dt_initial) + psi = self._get_avg_porosity(self._topo, solphi) + + self._solphi = solphi # Save solution of phi + self._psi_nm1 = psi # Average porosity value of last time step + + def _reinitialize_namespace(self, topo): + self._ns = None # Clear old namespace + self._ns = function.Namespace() + self._ns.x = self._geom + self._ns.ubasis = topo.basis('h-std', degree=self._degree_u).vector(topo.ndims) + self._ns.phibasis = topo.basis('h-std', degree=self._degree_phi) + self._ns.coarsephibasis = self._topo_coarse.basis('std', degree=self._degree_phi) + + self._ns.lam = (4 / self._nelems) / (2 ** self._ref_level) # Diffuse interface width + self._ns.gam = 0.05 + self._ns.kt = 1.0 + self._ns.eqconc = 0.5 # Equilibrium concentration + self._ns.kg = 0.0 # Conductivity of grain material + self._ns.ks = 1.0 # Conductivity of sand material + self._ns.reacrate = 'kt (?conc / eqconc)^2 - 1' # Constructed reaction rate based on macro temperature + self._ns.u = 'ubasis_ni ?solu_n' # Weights for which cell problem is solved for + self._ns.du_ij = 'u_i,j' # Gradient of weights field + self._ns.phi = 'phibasis_n ?solphi_n' # Phase field + self._ns.coarsephi = 'coarsephibasis_n ?coarsesolphi_n' # Phase field on original coarse topology + self._ns.ddwpdphi = '16 phi (1 - phi) (1 - 2 phi)' # gradient of double-well potential + self._ns.dphidt = 'phibasis_n (?solphi_n - ?solphinm1_n) / ?dt' # Implicit time evolution of phase field + + self._ucons = np.zeros(len(self._ns.ubasis), dtype=bool) + self._ucons[-1] = True # constrain u to zero at a point + + @staticmethod + def _analytical_phasefield(x, y, r, lam): + return 1. / (1. + function.exp(-4. / lam * (function.sqrt(x ** 2 + y ** 2) - r + 0.001))) + + @staticmethod + def _get_analytical_phasefield(topo, ns, degree_phi, lam, r): + phi_ini = MicroSimulation._analytical_phasefield(ns.x[0], ns.x[1], r, lam) + sqrphi = topo.integral((ns.phi - phi_ini) ** 2, degree=degree_phi * 2) + solphi = solver.optimize('solphi', sqrphi, droptol=1E-12) + + return solphi + + # def output(self): + # bezier = self._topo.sample('bezier', 2) + # x, u, phi = bezier.eval(['x_i', 'u_i', 'phi'] @ self._ns, solu=self._solu, solphi=self._solphi) + # with treelog.add(treelog.DataLog()): + # export.vtk("micro-heat", bezier.tri, x, T=u, phi=phi) + + def get_state(self): + return [self._solphi.copy(), deepcopy(self._topo)] + + def set_state(self, state): + self._solphi = state[0] + self._topo = state[1] + self._reinitialize_namespace(self._topo) # The namespace also needs to reloaded to its earlier state + + def _refine_mesh(self, topo_nm1, solphi_nm1): + """ + At the time of the calling of this function a predicted solution exists in ns.phi + """ + if not self._initial_condition_is_set: + solphi = solphi_nm1 + # ----- Refine the coarse mesh according to the projected solution to get a predicted refined topology ---- + topo = self._topo_coarse + for level in range(self._ref_level): + # print("refinement level = {}".format(level)) + smpl = topo.sample('uniform', 5) + ielem, criterion = smpl.eval([topo.f_index, abs(self._ns.coarsephi - .5) < .4], coarsesolphi=solphi) + + # Refine the elements for which at least one point tests true. + topo = topo.refined_by(np.unique(ielem[criterion])) + + self._reinitialize_namespace(topo) + + phi_ini = MicroSimulation._analytical_phasefield(self._ns.x[0], self._ns.x[1], self._r_initial, + self._ns.lam) + sqrphi = topo.integral((self._ns.coarsephi - phi_ini) ** 2, degree=self._degree_phi * 2) + solphi = solver.optimize('coarsesolphi', sqrphi, droptol=1E-12) + # ---------------------------------------------------------------------------------------------------- + else: + # ----- Refine the coarse mesh according to the projected solution to get a predicted refined topology ---- + topo = self._topo_coarse + for level in range(self._ref_level): + # print("refinement level = {}".format(level)) + topo_union1 = topo_nm1 & topo + smpl = topo_union1.sample('uniform', 5) + ielem, criterion = smpl.eval([topo.f_index, abs(self._ns.phi - .5) < .4], solphi=solphi_nm1) + + # Refine the elements for which at least one point tests true. + topo = topo.refined_by(np.unique(ielem[criterion])) + # ---------------------------------------------------------------------------------------------------- + + # Create a new projection mesh which is the union of the previous refined mesh and the predicted mesh + topo_union = topo_nm1 & topo + + # ----- Project the solution of the last time step on the projection mesh ----- + self._ns.projectedphi = function.dotarg('projectedsolphi', topo.basis('h-std', degree=self._degree_phi)) + sqrphi = topo_union.integral((self._ns.projectedphi - self._ns.phi) ** 2, degree=self._degree_phi * 2) + solphi = solver.optimize('projectedsolphi', sqrphi, droptol=1E-12, arguments=dict(solphi=solphi_nm1)) + + # Clip values of phase field to be in [0 1] to avoid overshoots and undershoots due to projection + solphi = np.clip(solphi, 0, 1) + + return topo, solphi + + def _solve_allen_cahn(self, topo, phi_coeffs_nm1, concentration, dt): + """ + Solving the Allen-Cahn equation using a Newton solver. + Returns porosity of the micro domain. + """ + self._first_iter_done = True + resphi = topo.integral('(lam^2 phibasis_n dphidt + gam phibasis_n ddwpdphi + gam lam^2 phibasis_n,i phi_,i + ' + '4 lam reacrate phibasis_n phi (1 - phi)) d:x' @ self._ns, degree=self._degree_phi * 2) + + args = dict(solphinm1=phi_coeffs_nm1, dt=dt, conc=concentration) + phi_coeffs = solver.newton('solphi', resphi, lhs0=phi_coeffs_nm1, arguments=args).solve(tol=1E-12) + + return phi_coeffs + + def _get_avg_porosity(self, topo, phi_coeffs): + psi = topo.integral('phi d:x' @ self._ns, degree=self._degree_phi * 2).eval(solphi=phi_coeffs) + + return psi + + def _solve_heat_cell_problem(self, topo, phi_coeffs): + """ + Solving the P1 homogenized heat equation + Returns upscaled conductivity for the micro domain + """ + res = topo.integral('((phi ks + (1 - phi) kg) u_i,j ubasis_ni,j - ' + '(ks - kg) phi_,j $_ij ubasis_ni) d:x' @ self._ns, degree=self._degree_u * 2) + + args = dict(solphi=phi_coeffs) + u_coeffs = solver.solve_linear('solu', res, constrain=self._ucons, arguments=args) + + return u_coeffs + + def _get_eff_conductivity(self, topo, u_coeffs, phi_coeffs): + b = topo.integral( + self._ns.eval_ij('(phi ks + (1 - phi) kg) ($_ij + du_ij) d:x'), + degree=self._degree_u * + 2).eval( + solu=u_coeffs, + solphi=phi_coeffs) + + return b.export("dense") + + def solve(self, macro_data, dt): + if self._psi_nm1 < 0.95: + topo, solphi = self._refine_mesh(self._topo, self._solphi) + self._reinitialize_namespace(topo) + + assert ((solphi >= 0.0) & (solphi <= 1.0)).all() + + solphi = self._solve_allen_cahn(topo, solphi, macro_data["concentration"], dt) + psi = self._get_avg_porosity(topo, solphi) + + solu = self._solve_heat_cell_problem(topo, solphi) + k = self._get_eff_conductivity(topo, solu, solphi) + + # Save state variables + self._topo = topo + self._solphi = solphi + self._solu = solu + self._psi_nm1 = psi + self._k_nm1 = k + else: + # Micro simulation has reached max porosity limit and hence is not solved + k = self._k_nm1 + psi = self._psi_nm1 + + output_data = dict() + output_data["k_00"] = k[0][0] + output_data["k_01"] = k[0][1] + output_data["k_10"] = k[1][0] + output_data["k_11"] = k[1][1] + output_data["porosity"] = psi + output_data["grain_size"] = math.sqrt((1 - psi) / math.pi) + + return output_data + + +def main(): + micro_problem = MicroSimulation() + dt = 1e-3 + micro_problem.initialize() + concentrations = [0.5, 0.4] + t = 0.0 + n = 0 + concentration = dict() + + for conc in concentrations: + concentration["concentration"] = conc + + micro_sim_output = micro_problem.solve(concentration, dt) + + # micro_problem.output() + t += dt + n += 1 + print(micro_sim_output) + + +if __name__ == "__main__": + cli.run(main) diff --git a/two-scale-heat-conduction/micro-nutils/run-micro-problems.py b/two-scale-heat-conduction/micro-nutils/run-micro-problems.py new file mode 100644 index 000000000..42900d3cf --- /dev/null +++ b/two-scale-heat-conduction/micro-nutils/run-micro-problems.py @@ -0,0 +1,11 @@ +""" +Script to run the Micro Manager +""" + +from micro_manager import MicroManager + +manager = MicroManager("micro-manager-config.json") + +manager.initialize() + +manager.solve() diff --git a/two-scale-heat-conduction/micro-nutils/run.sh b/two-scale-heat-conduction/micro-nutils/run.sh new file mode 100755 index 000000000..bff79b7f4 --- /dev/null +++ b/two-scale-heat-conduction/micro-nutils/run.sh @@ -0,0 +1,24 @@ +#!/bin/sh +set -e -u + +usage() { echo "Usage: cmd [-s] [-p n]" 1>&2; exit 1; } + +# Check if no input argument was provided +if [ -z "$*" ] ; then + echo "No input argument provided. Micro Manager is launched in serial" + python3 run-micro-problems.py +fi + +while getopts ":sp" opt; do + case ${opt} in + s) + python3 run-micro-problems.py + ;; + p) + mpiexec -n "$2" python3 run-micro-problems.py + ;; + *) + usage + ;; + esac +done diff --git a/two-scale-heat-conduction/precice-config.xml b/two-scale-heat-conduction/precice-config.xml new file mode 100644 index 000000000..df35681e2 --- /dev/null +++ b/two-scale-heat-conduction/precice-config.xml @@ -0,0 +1,104 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +