tasks/task08.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CAD integration\n",
    "\n",
    "In this task, we will learn how to run CAD-based simulations.\n",
    "\n",
    "Our example case will be a 3D ITER-like monoblock made of three different materials (tungsten, cucrzr, and copper).\n",
    "\n",
    "For this example, the mesh was created with [SALOME](https://www.salome-platform.org/) but other meshing software can be used like [GMSH](https://gmsh.info/).\n",
    "\n",
    "Using [meshio](https://github.com/nschloe/meshio), the mesh file (here a .med file) can be converted to .xdmf format, readable by FESTIM.\n",
    "\n",
    "We wrote the following function `convert_med_to_xdmf` to convert .med files to .xdmf."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import meshio\n",
    "\n",
    "\n",
    "def convert_med_to_xdmf(\n",
    "    med_file,\n",
    "    cell_file=\"mesh_domains.xdmf\",\n",
    "    facet_file=\"mesh_boundaries.xdmf\",\n",
    "    cell_type=\"tetra\",\n",
    "    facet_type=\"triangle\",\n",
    "):\n",
    "    \"\"\"Converts a MED mesh to XDMF\n",
    "    Args:\n",
    "        med_file (str): the name of the MED file\n",
    "        cell_file (str, optional): the name of the file containing the\n",
    "            volume markers. Defaults to \"mesh_domains.xdmf\".\n",
    "        facet_file (str, optional): the name of the file containing the\n",
    "            surface markers.. Defaults to \"mesh_boundaries.xdmf\".\n",
    "        cell_type (str, optional): The topology of the cells. Defaults to \"tetra\".\n",
    "        facet_type (str, optional): The topology of the facets. Defaults to \"triangle\".\n",
    "    Returns:\n",
    "        dict, dict: the correspondance dict, the cell types\n",
    "    \"\"\"\n",
    "    msh = meshio.read(med_file)\n",
    "\n",
    "    correspondance_dict = msh.cell_tags\n",
    "\n",
    "    cell_data_types = msh.cell_data_dict[\"cell_tags\"].keys()\n",
    "\n",
    "    for mesh_block in msh.cells:\n",
    "        if mesh_block.type == cell_type:\n",
    "\n",
    "            meshio.write_points_cells(\n",
    "                cell_file,\n",
    "                msh.points,\n",
    "                [mesh_block],\n",
    "                cell_data={\"f\": [-1 * msh.cell_data_dict[\"cell_tags\"][cell_type]]},\n",
    "            )\n",
    "        elif mesh_block.type == facet_type:\n",
    "            meshio.write_points_cells(\n",
    "                facet_file,\n",
    "                msh.points,\n",
    "                [mesh_block],\n",
    "                cell_data={\"f\": [-1 * msh.cell_data_dict[\"cell_tags\"][facet_type]]},\n",
    "            )\n",
    "\n",
    "    return correspondance_dict, cell_data_types"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The subdomains (volumes and surfaces) in the .med file were given tags.\n",
    "\n",
    "This is very important to be able to assign materials properties and boundary conditions to these subdomains.\n",
    "\n",
    "The correspondance between tags and subdomains is returned by `convert_med_to_xdmf`.\n",
    "Here for example, the volume corresponding to `tungsten` is the tag `6`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{-6: ['tungsten'], -7: ['cu'], -8: ['cucrzr'], -9: ['top_surface'], -10: ['cooling_surface'], -11: ['poloidal_gap'], -12: ['toroidal_gap'], -13: ['bottom']}\n"
     ]
    }
   ],
   "source": [
    "correspondance_dict, cell_data_types = convert_med_to_xdmf(\"task08/mesh.med\", cell_file=\"task08/mesh_domains.xdmf\", facet_file=\"task08/mesh_boundaries.xdmf\")\n",
    "\n",
    "print(correspondance_dict)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The converted .xdmf files can then be imported in FESTIM using the `MeshFromXDMF` class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Succesfully load mesh with 106966 cells\n"
     ]
    }
   ],
   "source": [
    "import festim as F\n",
    "\n",
    "my_model = F.Simulation()\n",
    "\n",
    "my_model.mesh = F.MeshFromXDMF(volume_file=\"task08/mesh_domains.xdmf\", boundary_file=\"task08/mesh_boundaries.xdmf\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the tags provided by `correspondance_dict`, we can create materials and assign them to the simulation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "tungsten = F.Material(\n",
    "    id=6,\n",
    "    D_0=4.1e-7,\n",
    "    E_D=0.39,\n",
    "    S_0=1.87e24,\n",
    "    E_S=1.04,\n",
    "    thermal_cond=100,\n",
    ")\n",
    "\n",
    "copper = F.Material(\n",
    "    id=7,\n",
    "    D_0=6.6e-7,\n",
    "    E_D=0.387,\n",
    "    S_0=3.14e24,\n",
    "    E_S=0.572,\n",
    "    thermal_cond=350,\n",
    ")\n",
    "\n",
    "cucrzr = F.Material(\n",
    "    id=8,\n",
    "    D_0=3.92e-7,\n",
    "    E_D=0.418,\n",
    "    S_0=4.28e23,\n",
    "    E_S=0.387,\n",
    "    thermal_cond=350\n",
    ")\n",
    "\n",
    "my_model.materials = [tungsten, copper, cucrzr]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarily, the surface tags are used to create boundary conditions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "heat_flux_top = F.FluxBC(surfaces=9, value=10e6, field=\"T\")\n",
    "convective_flux_coolant = F.ConvectiveFlux(surfaces=10, h_coeff=7e04, T_ext=323)\n",
    "implantation_flux_top = F.ImplantationDirichlet(surfaces=9, phi=1.61e22, R_p=9.52e-10, D_0=tungsten.D_0, E_D=tungsten.E_D)\n",
    "recombination_flux = F.DirichletBC(surfaces=[10, 11, 12, 13], value=0, field=0)\n",
    "\n",
    "my_model.boundary_conditions = [heat_flux_top, convective_flux_coolant, implantation_flux_top, recombination_flux]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will solve a steady state heat transfer problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_model.T = F.HeatTransferProblem(transient=False)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we add the settings (with `chemical_pot=True` to account for conservation of chemical potential) and we export the temperature and mobile concentration to XDMF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Defining variational problem heat transfers\n",
      "Solving stationary heat equation\n",
      "Calling FFC just-in-time (JIT) compiler, this may take some time.\n",
      "Calling FFC just-in-time (JIT) compiler, this may take some time.\n",
      "Calling FFC just-in-time (JIT) compiler, this may take some time.\n",
      "Defining initial values\n",
      "Defining variational problem\n",
      "Defining source terms\n",
      "Defining boundary conditions\n",
      "Solving steady state problem...\n",
      "Calling FFC just-in-time (JIT) compiler, this may take some time.\n",
      "  WARNING: The number of integration points for each cell will be: 125\n",
      "           Consider using the option 'quadrature_degree' to reduce the number of points\n",
      "  WARNING: The number of integration points for each cell will be: 125\n",
      "           Consider using the option 'quadrature_degree' to reduce the number of points\n",
      "  WARNING: The number of integration points for each cell will be: 125\n",
      "           Consider using the option 'quadrature_degree' to reduce the number of points\n",
      "Calling FFC just-in-time (JIT) compiler, this may take some time.\n",
      "  WARNING: The number of integration points for each cell will be: 125\n",
      "           Consider using the option 'quadrature_degree' to reduce the number of points\n",
      "  WARNING: The number of integration points for each cell will be: 125\n",
      "           Consider using the option 'quadrature_degree' to reduce the number of points\n",
      "  WARNING: The number of integration points for each cell will be: 125\n",
      "           Consider using the option 'quadrature_degree' to reduce the number of points\n",
      "Solved problem in 4.10 s\n"
     ]
    }
   ],
   "source": [
    "my_model.settings = F.Settings(\n",
    "    absolute_tolerance=1e10,\n",
    "    relative_tolerance=1e-10,\n",
    "    transient=False,\n",
    "    chemical_pot=True\n",
    ")\n",
    "\n",
    "results_folder = \"task08\"\n",
    "my_model.exports = [F.XDMFExport(\"T\", folder=results_folder), F.XDMFExport(\"solute\", folder=results_folder)]\n",
    "\n",
    "my_model.initialise()\n",
    "my_model.run()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Post processing"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The recommended way of plotting 3D fields is to export them to XDMF using `XDMFExport`, download the .xdmf file and open it in [Paraview](https://www.paraview.org/).\n",
    "\n",
    "This is what the fields look like in Paraview:"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img width=\"536\" alt=\"image\" src=\"https://github.com/pyvista/pyvista/assets/40028739/917a2d4f-96c0-4ad0-811f-7844e113c4aa\">\n",
    "<img width=\"520\" alt=\"image\" src=\"https://github.com/pyvista/pyvista/assets/40028739/324532f1-d5e9-42a5-8be4-8056da8c5f7c\">"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.12 ('base')",
   "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.8"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "d4d1e4263499bec80672ea0156c357c1ee493ec2b1c70f0acce89fc37c4a6abe"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}