Skip to content

Add plastic distortion tensor as output #2

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
tomswinburne wants to merge 2 commits into
base: main
Choose a base branch
from
Open

Add plastic distortion tensor as output #2

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
e7535ed
Add plastic distortion tensor output (Eptot, Wptot, Pdist)
tomswinburne Jun 5, 2026
eec1f94
Add doc on interpreting plastic distortion tensor output
tomswinburne Jun 5, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
65 changes: 65 additions & 0 deletions docs/plastic_distortion_output.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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.
45 changes: 43 additions & 2 deletions src/driver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand All @@ -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;
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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); }
Expand Down Expand Up @@ -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();
Expand Down Expand Up @@ -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");
Expand Down
6 changes: 5 additions & 1 deletion src/driver.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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;
Expand All @@ -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;
}
Expand Down