postProcess/getData-elastic-scalar2D.c
Simulation Data Extraction and Processing
This program extracts and processes data from fluid dynamics simulation snapshots, specifically designed for viscoelastic fluid simulations with conformation tensor analysis. It computes important derived quantities including deformation rate tensor components, velocity magnitude, and conformation tensor trace.
The program interpolates these quantities onto a regular grid and outputs the results for further analysis or visualization.
Physics Background
This code handles viscoelastic fluid simulation data where the fluid stress tensor includes both a viscous component (proportional to the deformation rate) and an elastic component (represented by the conformation tensor). The trace of the conformation tensor provides a measure of polymer stretching in the fluid.
#include "utils.h"
#include "output.h"
scalar f[]; // Volume fraction field
vector u[]; // Velocity field
scalar A11[], A12[], A22[]; // Conformation tensor components
scalar conform_qq[]; // Additional conformation tensor component
char filename[80]; // Input file name
int nx, ny, len; // Grid dimensions and field count
double xmin, ymin, xmax, ymax, Deltax, Deltay; // Domain boundaries and grid spacing
Derived Fields
- D2c: Log10 of squared deformation rate tensor weighted by volume fraction
- vel: Magnitude of velocity
- trA: Log10 of excess trace of conformation tensor
scalar D2c[], vel[], trA[];
scalar * list = NULL; // List to store output fields
Main Function
Processes simulation data and computes derived quantities
- Arguments:
- arguments[1]: Input filename
- arguments[2-5]: Domain boundaries (xmin, ymin, xmax, ymax)
- arguments[6]: Number of grid points in y-direction (ny)
- Returns:
- 0 on successful execution
int main(int a, char const *arguments[])
{
sprintf(filename, "%s", arguments[1]);
xmin = atof(arguments[2]); ymin = atof(arguments[3]);
xmax = atof(arguments[4]); ymax = atof(arguments[5]);
ny = atoi(arguments[6]);
// Initialize list of fields to output
list = list_add(list, D2c);
list = list_add(list, vel);
list = list_add(list, trA);
Data Processing Workflow
- Restore simulation state from snapshot file
- Compute derived quantities at each grid point
- Interpolate fields onto regular grid
- Output data to file
restore(file = filename);
Field Computation
For each cell, compute: - Components of the deformation rate tensor D - Squared magnitude of D weighted by volume fraction - Velocity magnitude - Excess trace of the conformation tensor
foreach() {
// Compute deformation rate tensor components
double D11 = (u.y[0, 1] - u.y[0, -1]) / (2 * Delta);
double D22 = (u.y[] / y);
double D33 = (u.x[1, 0] - u.x[-1, 0]) / (2 * Delta);
double D13 = 0.5 * ((u.y[1, 0] - u.y[-1, 0] + u.x[0, 1] - u.x[0, -1]) /
(2 * Delta));
double D2 = (sq(D11) + sq(D22) + sq(D33) + 2.0 * sq(D13));
D2c[] = f[] * D2;
// Take log10 of D2c for better visualization
if (D2c[] > 0.) {
D2c[] = log(D2c[]) / log(10);
} else {
D2c[] = -10;
}
// Compute velocity magnitude
vel[] = sqrt(sq(u.x[]) + sq(u.y[]));
// Compute excess trace of conformation tensor
trA[] = (A11[] + A22[] + conform_qq[]) / 3.0 - 1.0;
// Take log10 of trA for better visualization
if (trA[] > 0.) {
trA[] = log(trA[]) / log(10);
} else {
trA[] = -10;
}
}
Grid Interpolation and Output
- Calculate grid spacing based on domain size and ny
- Allocate memory for interpolated field values
- Interpolate field values onto regular grid
- Output grid coordinates and field values
FILE * fp = ferr;
Deltay = (double)((ymax - ymin) / (ny));
nx = (int)((xmax - xmin) / Deltay);
Deltax = (double)((xmax - xmin) / (nx));
len = list_len(list);
// Allocate memory for field values
double ** field = (double **) matrix_new(nx, ny + 1, len * sizeof(double));
// Interpolate field values onto regular grid
for (int i = 0; i < nx; i++) {
double x = Deltax * (i + 1./2) + xmin;
for (int j = 0; j < ny; j++) {
double y = Deltay * (j + 1./2) + ymin;
int k = 0;
for (scalar s in list) {
field[i][len * j + k++] = interpolate(s, x, y);
}
}
}
// Output grid coordinates and field values
for (int i = 0; i < nx; i++) {
double x = Deltax * (i + 1./2) + xmin;
for (int j = 0; j < ny; j++) {
double y = Deltay * (j + 1./2) + ymin;
fprintf(fp, "%g %g", x, y);
int k = 0;
for (scalar s in list) {
fprintf(fp, " %g", field[i][len * j + k++]);
}
fputc('\n', fp);
}
}
// Clean up
fflush(fp);
fclose(fp);
matrix_free(field);
return 0;
}