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. 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; }