{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cpfe-intro",
   "metadata": {},
   "source": [
    "# CPFE simulation\n",
    "\n",
    "Run a MOOSE/PUMA CPFE simulation on a far-field Voronoi reconstruction with a NEML2 v3 AOTI-compiled crystal-plasticity model, using `CPFESimulation`. Given a mesh and per-grain orientations (and an optional FF residual-strain field), it bakes the material parameters into the model, runs `neml2-compile`, and launches `puma-opt`.\n",
    "\n",
    "**Needs:** MOOSE `puma-opt` and `neml2-compile` (AOTI); CUDA recommended. See [installation](../install.rst).\n",
    "\n",
    "<div class=\"admonition warning\">\n",
    "<p class=\"admonition-title\">Tier B: needs MOOSE/PUMA</p>\n",
    "<p>This notebook needs <code>puma-opt</code> and <code>neml2-compile</code> (a GPU is recommended), which is not available in CI or Google Colab, so it is <strong>not executed when these docs are built</strong>. It ships as an annotated code listing with <strong>no stored output</strong> \u2014 nothing on this page is a hand-written or fabricated result. Run it in a configured environment (see the <a href=\"../install.rst\">installation guide</a>) to produce the figures, and run it from the <strong>repository root</strong> so the <code>mwe_data/</code> paths resolve.</p>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cpfe-inputs",
   "metadata": {},
   "source": "## Inputs\n\nA hex (`.e`) or tet (`.msh`) mesh plus per-grain orientations; CPFE runs on either. The repo ships `mwe_data/cpfe_ff/` (`reconstruction.msh` + `orientations.dat`, 10 grains). This sample is a GMSH tet mesh, but the recommended FF route is a SCULPT hex mesh (see [meshing](meshing.ipynb)). Orientations must be neml2 MRP; convert the FF Euler `orientations.dat` (degrees) with `orientation_helper.euler_to_mrp`."
  },
  {
   "cell_type": "markdown",
   "id": "cpfe-recipe-md",
   "metadata": {},
   "source": "## Recipe"
  },
  {
   "cell_type": "code",
   "id": "cpfe-recipe-code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "import os, numpy as np, torch, meshio\nfrom graintrace.run_cpfe_simulation import CPFESimulation\nfrom graintrace import orientation_helper as oh\n\nff = \"mwe_data/cpfe_ff\"; out = \"cpfe_out\"; os.makedirs(out, exist_ok=True)\neuler = np.loadtxt(ff + \"/orientations.dat\")                # Euler-bunge, degrees\nmrp = oh.euler_to_mrp(torch.tensor(euler, dtype=torch.float64), \"bunge\", \"degrees\")\nnp.savetxt(out + \"/orientations_MRP.dat\", mrp.numpy(), fmt=\"%.12g\")\n\nm = meshio.read(ff + \"/reconstruction.msh\"); P = m.points\nbbox = [float(P[:, 0].min()), float(P[:, 0].max()),\n        float(P[:, 1].min()), float(P[:, 1].max()),\n        float(P[:, 2].min()), float(P[:, 2].max())]\norder = \"SECOND\" if any(\"10\" in c.type for c in m.cells) else \"FIRST\"\n\nsim = CPFESimulation(\n    mesh_file=ff + \"/reconstruction.msh\", save_simulation_folder=out,\n    moose_run_file=\"external/puma/puma-opt\",   # your built puma-opt\n    element_order=order, eeres_file=None, ori_file=out + \"/orientations_MRP.dat\",\n    dim=3, use_ff_initial_field=True,           # eeres_file=None -> 12-col zero ee\n)\nsim.set_parameters(\"material\", slip_constant_strength=130.0,\n    voce_hardening_initial_slope=1556.09, voce_hardening_saturation=100.0,\n    power_slip_n=20, power_slip_g0=1e-4, elastic_E=209016.0, elastic_nu=0.307,\n    elastic_G=60355.0, burger_scale=2.22)\nsim.set_parameters(\"simulation_parameters\", device=\"cuda:0\", device_batch=20000,\n    dt=0.5, total_time=2.0, initialize_time=1.0, sync_times=\"2.0\",\n    grid_transfer=\"final\", exodus_output=\"sync\",\n    distributed_mesh=False)  # True -> pre-split (--split-mesh) + --use-split; large meshes, ncore>=2\n\ndisp = 0.002 * (bbox[5] - bbox[4])\nsim.set_parameters(\"boundary\", bounding_box=bbox, bc={\n    \"x\": {\"negative\": \"stress_free\", \"positive\": \"stress_free\"},\n    \"y\": {\"negative\": \"stress_free\", \"positive\": \"stress_free\"},\n    \"z\": {\"negative\": 0, \"positive\": disp}})\ngrid_bb = list(bbox)\nfor i in range(0, 6, 2): grid_bb[i] += 1e-4\nfor i in range(1, 6, 2): grid_bb[i] -= 1e-4\nsim.set_parameters(\"grid_properties\", number_of_elements=[10, 10, 10],\n                   bounding_box=grid_bb)\n\nsim.run(ncore=4)   # ncore == mpiexec -n; also spreads a device list over ranks"
  },
  {
   "cell_type": "markdown",
   "id": "cpfe-key-params",
   "metadata": {},
   "source": "## Key parameters\n\n- `device`: `\"cpu\"`, `\"cuda:0\"`, or a space-separated list `\"cuda:0 cuda:1\"` (multi-GPU over MPI ranks).\n- `device_batch`: per-device NEML2 chunk (quad points per call); a finite value caps GPU memory (0 means the whole batch, risking OOM on large meshes).\n- `initialize_time`: the load ramps from `initialize_time` to `total_time`; `sync_times` are the grid-output times.\n- `use_ff_initial_field=True` with a real `eeres_file` applies an FF residual strain (12 columns: x, y, z plus 9); `eeres_file=None` writes a 12-column zero ee.\n- `grid_transfer` and `exodus_output` default to the cheap settings (transfer only at the last step, Exodus only at `sync_times`). Use `\"per_step\"` for every-step grid REI output or Exodus writes; with `grid_transfer` other than `\"per_step\"`, regenerate grid CSVs offline with `GridResampler`.\n- `distributed_mesh`: `True` pre-splits the mesh (`--split-mesh ncore`) and runs `--use-split` so each rank reads only its partition \u2014 the memory lever for large meshes (replicated mesh holds the full mesh on every rank). Pre-split only (no in-situ option); requires `ncore >= 2` and a `puma-opt` build with the `EqualValueBoundaryConstraint` distributed-mesh fix. Outputs are unchanged (single `sim_output.e`, complete `mesh_out/`/`grid_out/`). Default `False`.\n- `solver_route`: `\"timestep_optimized\"` (default) uses a direct LU solve (`superlu_dist`) with preconditioner reuse \u2014 robust and iteration-cheap, but memory-bound on large meshes; `\"hpc_memory\"` uses an iterative `fgmres` + GAMG algebraic-multigrid solve \u2014 low, scalable memory for large HPC meshes and the recommended partner of `distributed_mesh=True`, at the cost of more/looser linear iterations. Both keep `residual_and_jacobian_together = false` (nodal-constraint loading BC; MOOSE issue 33531).\n- AOTI: material parameters are baked into the model `.i` and `neml2-compile`d on `run()`; `recompile=True` (default) rebuilds when parameters change.\n\nSee [configuration](../configuration.rst) for the `bc` dict and inset `grid_properties` bounding box."
  },
  {
   "cell_type": "markdown",
   "id": "89c1e421",
   "source": "## Running on HPC\n\nFor a cluster run, switch the launcher to the scheduler and use the memory-lean solver; add the distributed mesh only once the replicated mesh no longer fits:\n\n```python\nsim.set_parameters(\"simulation_parameters\",\n    solver_route=\"hpc_memory\",   # iterative fgmres + GAMG (low, scalable memory)\n    launcher=\"srun\",             # Slurm/Cray; default \"mpiexec\" elsewhere\n    distributed_mesh=True,       # large meshes only (~1M+ elements); needs ncore >= 2\n    device=\"cuda:0 cuda:1 cuda:2 cuda:3\",  # GPU: a cuda list, one entry per GPU on the node\n    device_batch=20000)          # finite chunk caps per-GPU memory (0 risks OOM)\nsim.run(ncore=4)                 # GPU: ncore == number of GPUs in the device list\n```\n\n- **GPU node:** `device` is a space-separated cuda list and `ncore` equals the number of GPUs; fewer, larger ranks give better GPU utilisation for the NEML2-dominated solve. Keep `device_batch` finite to bound per-GPU memory.\n- **CPU node:** set `device=\"cpu\"` and `ncore` to the number of MPI ranks (`srun -n`); `distributed_mesh=True` is what keeps per-rank memory bounded as the mesh grows.\n- `solver_route=\"hpc_memory\"` and `distributed_mesh=True` are independent levers but pair naturally \u2014 GAMG avoids the LU factorisation blow-up, the distributed mesh avoids holding the full mesh per rank. Enable `distributed_mesh` only for meshes that OOM as replicated (needs the EVBC-fixed `puma-opt`); `hpc_memory` is safe at any size.",
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "id": "cpfe-outputs",
   "metadata": {},
   "source": "## Outputs (in `save_simulation_folder/simulation_out`)\n\n- `out.csv`: per-block time series.\n- `sim_output.e` / `sim_output_grid.e`: Exodus output.\n- `grid_out/*.csv`: per-grid fields (cauchy_stress, ee, nye_tensor, ori_rodrigues), for post-processing and REI."
  },
  {
   "cell_type": "markdown",
   "id": "cpfe-gotchas",
   "metadata": {},
   "source": "## Gotchas\n\n- v3 has no runtime `[NEML2] cli_args` or `[Schedulers]`; do not pass `scheduler_name`.\n- For neml2-dominated CPFE, use fewer MPI ranks (about one per GPU) for bigger per-rank batches.\n- Stiff first steps may make MOOSE cut `dt` and recover; that is normal, not a failure.\n- The environment, `LD_LIBRARY_PATH`, and `neml2_load_files` auto-derive from the `moose_run_file` repo layout."
  },
  {
   "cell_type": "markdown",
   "id": "cpfe-full-example",
   "metadata": {},
   "source": "## Full example\n\n**Full script:** [examples/demonstrate_cpfe.py](https://github.com/applied-material-modeling/graintrace/blob/main/examples/demonstrate_cpfe.py)"
  }
 ],
 "metadata": {
  "nbsphinx": {
   "execute": "never"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}