From e7535ed84b0749191d52771986e2a0379cc04b17 Mon Sep 17 00:00:00 2001 From: tomswinburne Date: Fri, 5 Jun 2026 13:40:42 -0400 Subject: [PATCH 1/2] Add plastic distortion tensor output (Eptot, Wptot, Pdist) Accumulate the symmetric plastic strain (Eptot) and plastic spin (Wptot) per step in the driver, expose them along with the full plastic distortion tensor (Pdist = Eptot + Wptot) as output fields, and persist Eptot/Wptot across restarts. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/driver.cpp | 45 +++++++++++++++++++++++++++++++++++++++++++-- src/driver.h | 6 +++++- 2 files changed, 48 insertions(+), 3 deletions(-) diff --git a/src/driver.cpp b/src/driver.cpp index b9b81a0..c8dde0d 100644 --- a/src/driver.cpp +++ b/src/driver.cpp @@ -25,6 +25,8 @@ ExaDiSApp::ExaDiSApp(int argc, char* argv[]) istep = 0; Etot = Mat33().zero(); + Eptot = Mat33().zero(); + Wptot = Mat33().zero(); stress = strain = pstrain = 0.0; tottime = 0.0; @@ -46,6 +48,8 @@ ExaDiSApp::ExaDiSApp() { istep = 0; Etot = Mat33().zero(); + Eptot = Mat33().zero(); + Wptot = Mat33().zero(); stress = strain = pstrain = 0.0; tottime = 0.0; dealloc = false; @@ -160,6 +164,10 @@ void ExaDiSApp::write_restart(std::string restartfile) fprintf(fp, "step %d\n", istep); fprintf(fp, "Etot %.17g %.17g %.17g %.17g %.17g %.17g %.17g %.17g %.17g\n", Etot.xx(), Etot.xy(), Etot.xz(), Etot.yx(), Etot.yy(), Etot.yz(), Etot.zx(), Etot.zy(), Etot.zz()); + fprintf(fp, "Eptot %.17g %.17g %.17g %.17g %.17g %.17g %.17g %.17g %.17g\n", + Eptot.xx(), Eptot.xy(), Eptot.xz(), Eptot.yx(), Eptot.yy(), Eptot.yz(), Eptot.zx(), Eptot.zy(), Eptot.zz()); + fprintf(fp, "Wptot %.17g %.17g %.17g %.17g %.17g %.17g %.17g %.17g %.17g\n", + Wptot.xx(), Wptot.xy(), Wptot.xz(), Wptot.yx(), Wptot.yy(), Wptot.yz(), Wptot.zx(), Wptot.zy(), Wptot.zz()); fprintf(fp, "stress %.17g\n", stress); fprintf(fp, "strain %.17g\n", strain); fprintf(fp, "pstrain %.17g\n", pstrain); @@ -264,6 +272,18 @@ void ExaDiSApp::read_restart(std::string restartfile) &Etot[1][0], &Etot[1][1], &Etot[1][2], &Etot[2][0], &Etot[2][1], &Etot[2][2]); } + else if (strncmp(line, "Eptot", 5) == 0) { + sscanf(line, "Eptot %lf %lf %lf %lf %lf %lf %lf %lf %lf\n", + &Eptot[0][0], &Eptot[0][1], &Eptot[0][2], + &Eptot[1][0], &Eptot[1][1], &Eptot[1][2], + &Eptot[2][0], &Eptot[2][1], &Eptot[2][2]); + } + else if (strncmp(line, "Wptot", 5) == 0) { + sscanf(line, "Wptot %lf %lf %lf %lf %lf %lf %lf %lf %lf\n", + &Wptot[0][0], &Wptot[0][1], &Wptot[0][2], + &Wptot[1][0], &Wptot[1][1], &Wptot[1][2], + &Wptot[2][0], &Wptot[2][1], &Wptot[2][2]); + } else if (strncmp(line, "stress", 6) == 0) { sscanf(line, "stress %lf\n", &stress); } else if (strncmp(line, "strain", 6) == 0) { sscanf(line, "strain %lf\n", &strain); } else if (strncmp(line, "pstrain", 7) == 0) { sscanf(line, "pstrain %lf\n", &pstrain); } @@ -382,6 +402,11 @@ void ExaDiSApp::read_restart(std::string restartfile) *-------------------------------------------------------------------------*/ void ExaDiSApp::update_mechanics(Control& ctrl) { + // Accumulate the full plastic distortion tensor from the per-step + // increments (Eptot = symmetric strain, Wptot = antisymmetric spin) + Eptot += system->dEp; + Wptot += system->dWp; + if (ctrl.rotation) { // Counter-rotate stress wrt dislocation configuration double p1 = - system->dWp.zy(); @@ -479,9 +504,25 @@ void ExaDiSApp::output(Control& ctrl) if (header) fprintf(fp, " Rxx Rxy Rxz Ryx Ryy Ryz Rzx Rzy Rzz"); else fprintf(fp, "%e %e %e %e %e %e %e %e %e ", R.xx(), R.xy(), R.xz(), R.yx(), R.yy(), R.yz(), R.zx(), R.zy(), R.zz()); } else if (field == Prop::ALLSTRESS) { - if (header) fprintf(fp, " Sxx Syy Szz Sxy Sxz Syz"); - else fprintf(fp, "%e %e %e %e %e %e ", system->extstress.xx(), system->extstress.yy(), system->extstress.zz(), + if (header) fprintf(fp, " Sxx Syy Szz Sxy Sxz Syz"); + else fprintf(fp, "%e %e %e %e %e %e ", system->extstress.xx(), system->extstress.yy(), system->extstress.zz(), system->extstress.xy(), system->extstress.xz(), system->extstress.yz()); + } else if (field == Prop::EPTOT) { + if (header) fprintf(fp, " Epxx Epxy Epxz Epyx Epyy Epyz Epzx Epzy Epzz"); + else fprintf(fp, "%e %e %e %e %e %e %e %e %e ", Eptot.xx(), Eptot.xy(), Eptot.xz(), + Eptot.yx(), Eptot.yy(), Eptot.yz(), + Eptot.zx(), Eptot.zy(), Eptot.zz()); + } else if (field == Prop::WPTOT) { + if (header) fprintf(fp, " Wpxx Wpxy Wpxz Wpyx Wpyy Wpyz Wpzx Wpzy Wpzz"); + else fprintf(fp, "%e %e %e %e %e %e %e %e %e ", Wptot.xx(), Wptot.xy(), Wptot.xz(), + Wptot.yx(), Wptot.yy(), Wptot.yz(), + Wptot.zx(), Wptot.zy(), Wptot.zz()); + } else if (field == Prop::PDIST) { + Mat33 Bp = Eptot + Wptot; // full plastic distortion tensor + if (header) fprintf(fp, " Bpxx Bpxy Bpxz Bpyx Bpyy Bpyz Bpzx Bpzy Bpzz"); + else fprintf(fp, "%e %e %e %e %e %e %e %e %e ", Bp.xx(), Bp.xy(), Bp.xz(), + Bp.yx(), Bp.yy(), Bp.yz(), + Bp.zx(), Bp.zy(), Bp.zz()); } } fprintf(fp, "\n"); diff --git a/src/driver.h b/src/driver.h index e32a2f6..8ce5690 100644 --- a/src/driver.h +++ b/src/driver.h @@ -42,6 +42,7 @@ class ExaDiSApp { int istep; Mat33 Etot; + Mat33 Eptot, Wptot; // accumulated plastic strain and plastic spin tensors double stress, strain, pstrain; double tottime; Vec3 edir; @@ -68,7 +69,7 @@ class ExaDiSApp { static Stepper MAX_WALLTIME(double maxtime) { return Stepper(Stepper::MAX_WALLTIME, maxtime); } struct Prop { - enum fields {STEP, STRAIN, STRESS, DENSITY, NNODES, NSEGS, DT, TIME, WALLTIME, EDIR, RORIENT, ALLSTRESS}; + enum fields {STEP, STRAIN, STRESS, DENSITY, NNODES, NSEGS, DT, TIME, WALLTIME, EDIR, RORIENT, ALLSTRESS, EPTOT, WPTOT, PDIST}; static fields get_field(std::string& name) { if (name == "Step" || name == "step") return STEP; else if (name == "Strain" || name == "strain") return STRAIN; @@ -82,6 +83,9 @@ class ExaDiSApp { else if (name == "edir") return EDIR; else if (name == "Rorient") return RORIENT; else if (name == "Allstress" || name == "allstress") return ALLSTRESS; + else if (name == "Eptot" || name == "eptot") return EPTOT; + else if (name == "Wptot" || name == "wptot") return WPTOT; + else if (name == "Pdist" || name == "pdist") return PDIST; else ExaDiS_fatal("Unknown control property name = %s\n", name.c_str()); return STEP; } From eec1f940efd21c8881e8777bcaee1909e38178a0 Mon Sep 17 00:00:00 2001 From: tomswinburne Date: Fri, 5 Jun 2026 14:22:45 -0400 Subject: [PATCH 2/2] Add doc on interpreting plastic distortion tensor output Co-Authored-By: Claude Opus 4.8 (1M context) --- docs/plastic_distortion_output.md | 65 +++++++++++++++++++++++++++++++ 1 file changed, 65 insertions(+) create mode 100644 docs/plastic_distortion_output.md diff --git a/docs/plastic_distortion_output.md b/docs/plastic_distortion_output.md new file mode 100644 index 0000000..1c88200 --- /dev/null +++ b/docs/plastic_distortion_output.md @@ -0,0 +1,65 @@ +# Plastic distortion tensor output + +Adds three accumulated plastic-deformation output fields to the driver: +`Eptot`, `Wptot`, and `Pdist`. They are summed from the per-step increments +(`system->dEp`, `system->dWp`) in `ExaDiSApp::update_mechanics`, persisted in +restart files, and emitted as output columns when requested. + +## The fields + +| Field | Meaning | Symmetry | Columns | +|---------|----------------------------------|---------------|---------| +| `Eptot` | Accumulated plastic strain | symmetric | 9 (full 3×3) | +| `Wptot` | Accumulated plastic spin/rotation| antisymmetric | 9 (full 3×3) | +| `Pdist` | Full plastic distortion `βp` | general | 9 (full 3×3) | + +Relationship: + +``` +Pdist = Eptot + Wptot (βp = εp + ωp) +Eptot = ½ (βp + βpᵀ) symmetric part → plastic strain +Wptot = ½ (βp − βpᵀ) antisymmetric part → plastic rotation +``` + +## Column layout + +Each field is written as a full 3×3 tensor in **row-major** order: + +``` +xx xy xz yx yy yz zx zy zz +``` + +So the header for `Eptot` is `Epxx Epxy Epxz Epyx Epyy Epyz Epzx Epzy Epzz` +(`Wp…` for `Wptot`, `Bp…` for `Pdist`). + +## Enabling the output + +Add the property name to the driver's output property list (case-insensitive): + +| Property name | Field | +|----------------------|---------| +| `Eptot` / `eptot` | `Eptot` | +| `Wptot` / `wptot` | `Wptot` | +| `Pdist` / `pdist` | `Pdist` | + +## Interpreting values + +- **`Eptot` (plastic strain).** Symmetric. Diagonal = normal plastic strains + along x/y/z; off-diagonal = engineering shear (½γ) components. For a single + active slip system you'll see the resolved shear dominate one off-diagonal + pair. Trace ≈ 0 — dislocation glide is volume-preserving (plastic + incompressibility). +- **`Wptot` (plastic spin).** Antisymmetric (zero diagonal, `ij = −ji`). It is + the lattice rotation accumulated by plastic flow — relevant for texture + evolution and for the stress counter-rotation applied when + `ctrl.rotation` is on. +- **`Pdist` (`βp`).** The full displacement-gradient contribution from + plasticity. Use it when you need strain and rotation together (e.g. comparing + against a continuum `βp`), rather than reconstructing it from the two halves. + +## Restart + +`Eptot` and `Wptot` are written to and read back from the restart file +(`write_restart` / `read_restart`), so accumulation continues seamlessly across +restarts. `Pdist` is derived on the fly (`Eptot + Wptot`) and is not stored +separately.