Skip to content

semi-optimized for a 3d structure #1

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
wband wants to merge 1 commit into
base: main
Choose a base branch
from
Open

semi-optimized for a 3d structure #1

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
85 changes: 85 additions & 0 deletions grains_in_c++/3D_grain_structure/datread.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,85 @@
import numpy as np
import matplotlib.pyplot as plt
plt.set_cmap('Greys')
nx = 65
ny = 65
nz = 65
imin = np.fromfile("grain_vis.dat")
#imin = np.fromfile("30571_seg.dat", offset=0)
#imin = np.fromfile("ACnon_100.dat", offset=4)
#imin = np.fromfile("ACcon_non400.dat", offset=4)
#imin = np.fromfile("piecestest/ACcon_big.dat", offset=4)
#imin = np.fromfile("piecestest/in0002.dat", offset=4)
print(len(imin))
im = imin.reshape([nx,ny,nz],order='F')
#im = imin.reshape([150,173,90],order='F')
#im = imin.reshape([150,345,360],order='F')

print(im.shape)
print(im[0,0,0],im[-1,-1,-1])
print(np.amax(imin),np.amin(imin))

stem='CHds2_t0010k'
#stem = 'AC100'
# plot faces
vmin_all = 0
vmax_all = 1.0
fig = plt.figure(figsize=(3,3))
fig.patch.set_visible(False)
ax = fig.add_axes((0,0,1,1))
ax.set_axis_off()
axim = ax.imshow(im[0,:,:], origin='lower', vmin=vmin_all, vmax=vmax_all)
fig.savefig(stem+'x0.pdf',dpi=300)
fig.clear()
ax = fig.add_axes((0,0,1,1))
ax.set_axis_off()
axim = ax.imshow(im[:,0,:], origin='lower', vmin=vmin_all, vmax=vmax_all)
fig.savefig(stem+'y0.pdf',dpi=300)
fig.clear()
ax = fig.add_axes((0,0,1,1))
ax.set_axis_off()
axim = ax.imshow(im[:,:,0], origin='lower', vmin=vmin_all, vmax=vmax_all)
fig.savefig(stem+'z0.pdf',dpi=300)
fig.clear()

ax = fig.add_axes((0,0,1,1))
ax.set_axis_off()
axim = ax.imshow(im[-1,:,:], origin='lower', vmin=vmin_all, vmax=vmax_all)
fig.savefig(stem+'x1.pdf',dpi=300)
fig.clear()
ax = fig.add_axes((0,0,1,1))
ax.set_axis_off()
axim = ax.imshow(im[:,-1,:], origin='lower', vmin=vmin_all, vmax=vmax_all)
fig.savefig(stem+'y1.pdf',dpi=300)
fig.clear()
ax = fig.add_axes((0,0,1,1))
ax.set_axis_off()
axim = ax.imshow(im[:,:,-1], origin='lower', vmin=vmin_all, vmax=vmax_all)
fig.savefig(stem+'z1.pdf',dpi=300)
fig.clear()

ax = fig.add_axes((0,0,1,1))
ax.set_axis_off()
axim = ax.imshow(im[int(nx/2),:,:], origin='lower', vmin=vmin_all, vmax=vmax_all)
fig.savefig(stem+'x12.pdf',dpi=300)
fig.clear()
ax = fig.add_axes((0,0,1,1))
ax.set_axis_off()
axim = ax.imshow(im[int(nx/2)+1,:,:], origin='lower', vmin=vmin_all, vmax=vmax_all)
fig.savefig(stem+'x12_1.pdf',dpi=300)
fig.clear()
ax = fig.add_axes((0,0,1,1))
ax.set_axis_off()
axim = ax.imshow(im[int(nx/2)+2,:,:], origin='lower', vmin=vmin_all, vmax=vmax_all)
fig.savefig(stem+'x12_2.pdf',dpi=300)
fig.clear()
ax = fig.add_axes((0,0,1,1))
ax.set_axis_off()
axim = ax.imshow(im[:,int(ny/2),:], origin='lower', vmin=vmin_all, vmax=vmax_all)
fig.savefig(stem+'x13.pdf',dpi=300)
fig.clear()
ax = fig.add_axes((0,0,1,1))
ax.set_axis_off()
axim = ax.imshow(im[:,:,int(nz/2)], origin='lower', vmin=vmin_all, vmax=vmax_all)
fig.savefig(stem+'x23.pdf',dpi=300)
fig.clear()
137 changes: 77 additions & 60 deletions grains_in_c++/3D_grain_structure/grain_maker.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,17 +10,19 @@ using namespace std;

int main() {
// --- PARAMETERS ---
const int nx = 128, ny = 128, nz = 128;
const int nx = 97, ny = 97, nz = 97;
const int pad_x = nx + 2, pad_y = ny + 2,
pad_z = nz + 2; // Including buffer layers (x1-1 to x2+1)
const int np_initial = 100;
const int ssteps = 500;
const int extrax = 16;
const int innerx = nx - 2*extrax, innery = ny - 2*extrax, innerz = nz - 2*extrax;
const int np_initial = 200;
const int ssteps = 700;

const double dh = 1.0;
const double gamma = 1.5;
const double W = 1.0;
const double eps_sq = 1.0;
const double sLdt = 0.1;
const double sLdt = 0.08;

// Helper lambdas for multi-dimensional 1D indexing
// Layout: [z][y][x][p] to optimize cache locality during the spatial stencils
Expand All @@ -31,6 +33,14 @@ int main() {
auto idx4 = [pad_x, pad_y](int x, int y, int z, int p, int num_p) {
return ((z * pad_y * pad_x) + (y * pad_x) + x) * num_p + p;
};

auto idx3_small = [innerx, innery](int x, int y, int z) {
return (z * innery * innerx) + (y * innerx) + x;
};

auto idx4_small = [innerx, innery](int x, int y, int z, int p, int num_p) {
return ((z * innery * innerx) + (y * innerx) + x) * num_p + p;
};

// --- INITIAL VORONOI TESSELLATION ---
cout << "Initializing Voronoi tessellation..." << endl;
Expand All @@ -46,6 +56,11 @@ int main() {
centers[ip][1] = dist(gen) * length[1];
centers[ip][2] = dist(gen) * length[2];
}

// ip = 0 is now the big grain
centers[0][0] = 0.0;
centers[0][1] = 0.0;
centers[0][2] = 0.0;

for (int iz = 0; iz < pad_z; ++iz) {
for (int iy = 0; iy < pad_y; ++iy) {
Expand All @@ -69,25 +84,20 @@ int main() {
}
}
}

// --- KNOCK OUT GRAINS ---
cout << "Knocking out invalid grains..." << endl;
vector<int> index_to_featureID;
int npeff = 0;

for (int ip = 0; ip < np_initial; ++ip) {
for (int ip = 1; ip < np_initial; ++ip) {
int count_mask = 0;
for (int v : featureID) {
if (v == ip)
count_mask++;
}

double dist2_center = pow(centers[ip][0] - (nx / 2.0), 2) +
pow(centers[ip][1] - (ny / 2.0), 2) +
pow(centers[ip][2] - (nz / 2.0), 2);
double radius_limit = pow((nx / 2.0) * 0.85, 2);

if (count_mask < 300 || dist2_center > radius_limit) {
if (count_mask < 100) {
for (int &v : featureID) {
if (v == ip)
v = np_initial; // mapped to invalid grain
Expand All @@ -106,12 +116,17 @@ int main() {

// --- INITIALIZE ORDER PARAMETERS (PHI) ---
vector<double> Phi(pad_x * pad_y * pad_z * npeff, 0.0);
for (int ip = 0; ip < npeff; ++ip) {
for (int iz = 0; iz < pad_z; ++iz) {
for (int iy = 0; iy < pad_y; ++iy) {
for (int ix = 0; ix < pad_x; ++ix) {
if (featureID[idx3(ix, iy, iz)] == index_to_featureID[ip]) {
Phi[idx4(ix, iy, iz, ip, npeff)] = 1.0;
for (int iz = 0; iz < pad_z; ++iz) {
for (int iy = 0; iy < pad_y; ++iy) {
for (int ix = 0; ix < pad_x; ++ix) {
double dist = sqrt(pow(ix - (nx / 2.0), 2) +
pow(iy - (ny / 2.0), 2) +
pow(iz - (nz / 2.0), 2));
double distfn = 0.5*(1.0+tanh((innerx/2*1.2-dist)/2.0));
Phi[idx4(ix, iy, iz, 0, npeff)] = 1.0-distfn;
for (int ip = 0; ip < npeff; ++ip) {
if (featureID[idx3(ix, iy, iz)] == ip) {
Phi[idx4(ix, iy, iz, ip, npeff)] = distfn;
}
}
}
Expand All @@ -120,10 +135,32 @@ int main() {

featureID.clear(); // Free memory
index_to_featureID.clear();

// Applying Dirichlet BCs
for (int ip = 0; ip < npeff; ++ip) {
for (int iz = 0; iz < pad_z; ++iz) {
for (int iy = 0; iy < pad_y; ++iy) {
Phi[idx4(0, iy, iz, ip, npeff)] = 0.0;
Phi[idx4(nx + 1, iy, iz, ip, npeff)] = 0.0;
}
for (int ix = 0; ix < pad_x; ++ix) {
Phi[idx4(ix, 0, iz, ip, npeff)] = 0.0;
Phi[idx4(ix, ny + 1, iz, ip, npeff)] = 0.0;
}
}
for (int iy = 0; iy < pad_y; ++iy) {
for (int ix = 0; ix < pad_x; ++ix) {
Phi[idx4(ix, iy, 0, ip, npeff)] = 0.0;
Phi[idx4(ix, iy, nz + 1, ip, npeff)] = 0.0;
}
}
}

vector<double> Phi_new = Phi;
vector<double> smsq(pad_x * pad_y * pad_z, 0.0);



// --- TIME EVOLUTION ---
cout << "Starting time integration (" << ssteps << " steps)..." << endl;
for (int it = 0; it <= ssteps; ++it) {
Expand Down Expand Up @@ -167,32 +204,6 @@ int main() {
}
}

// Apply Periodic BCs in buffer layers
for (int ip = 0; ip < npeff; ++ip) {
for (int iz = 0; iz < pad_z; ++iz) {
for (int iy = 0; iy < pad_y; ++iy) {
Phi_new[idx4(0, iy, iz, ip, npeff)] =
Phi_new[idx4(nx, iy, iz, ip, npeff)];
Phi_new[idx4(nx + 1, iy, iz, ip, npeff)] =
Phi_new[idx4(1, iy, iz, ip, npeff)];
}
for (int ix = 0; ix < pad_x; ++ix) {
Phi_new[idx4(ix, 0, iz, ip, npeff)] =
Phi_new[idx4(ix, ny, iz, ip, npeff)];
Phi_new[idx4(ix, ny + 1, iz, ip, npeff)] =
Phi_new[idx4(ix, 1, iz, ip, npeff)];
}
}
for (int iy = 0; iy < pad_y; ++iy) {
for (int ix = 0; ix < pad_x; ++ix) {
Phi_new[idx4(ix, iy, 0, ip, npeff)] =
Phi_new[idx4(ix, iy, nz, ip, npeff)];
Phi_new[idx4(ix, iy, nz + 1, ip, npeff)] =
Phi_new[idx4(ix, iy, 1, ip, npeff)];
}
}
}

Phi = Phi_new;

// --- DYNAMIC GRAIN PRUNING ---
Expand Down Expand Up @@ -231,12 +242,15 @@ int main() {
Phi = Phi_shrink;
Phi_new = Phi;
npeff = npeff_fin;
} else if (it % 500 == 0) {
cout << "Step " << it << " completed." << endl;
}
}
if (it % 500 == 0) {
cout << "Step " << it << " completed." << endl;
}
}



// --- POST-PROCESSING ---
cout << "Normalizing and writing outputs..." << endl;
for (double &val : Phi) {
Expand Down Expand Up @@ -283,13 +297,18 @@ int main() {
// Extract interior grid for outputs (ignoring pad layers, mapping to exact
// dimensions make_props expects)

cout << "writing into " << innerx << "x" << innery << "x" << innerz << endl;

// 1. grain_vis.dat
vector<double> smsq_out(nx * ny * nz);
for (int iz = 1; iz <= nz; ++iz)
for (int iy = 1; iy <= ny; ++iy)
for (int ix = 1; ix <= nx; ++ix)
smsq_out[((iz - 1) * ny + (iy - 1)) * nx + (ix - 1)] =
smsq[idx3(ix, iy, iz)];
vector<double> smsq_out(innerx * innery * innerz);
for (int iz = 0; iz <= innerz-1; ++iz) {
for (int iy = 0; iy <= innery-1; ++iy) {
for (int ix = 0; ix <= innerx-1; ++ix) {
smsq_out[idx3_small(ix,iy,iz)] =
smsq[idx3(ix+extrax, iy+extrax, iz+extrax)];
}
}
}

ofstream outVis("grain_vis.dat", ios::binary);
if (outVis.is_open()) {
Expand All @@ -300,14 +319,12 @@ int main() {

// 2. grains.dat (This must perfectly match the 4D interior space for your
// other C++ script)
vector<double> Phi_out(nx * ny * nz * npeff);
for (int iz = 1; iz <= nz; ++iz) {
for (int iy = 1; iy <= ny; ++iy) {
for (int ix = 1; ix <= nx; ++ix) {
vector<double> Phi_out(innerx * innery * innerz * npeff);
for (int iz = 0; iz <= innerz-1; ++iz) {
for (int iy = 0; iy <= innery-1; ++iy) {
for (int ix = 0; ix <= innerx-1; ++ix) {
for (int ip = 0; ip < npeff; ++ip) {
int out_idx = (((ix - 1) * ny + (iy - 1)) * nz + (iz - 1)) * npeff +
ip; // matches idx4D in make_props
Phi_out[out_idx] = Phi[idx4(ix, iy, iz, ip, npeff)];
Phi_out[idx4_small(ix,iy,iz,ip, npeff)] = Phi[idx4(ix+extrax, iy+extrax, iz+extrax, ip, npeff)];
}
}
}
Expand All @@ -322,4 +339,4 @@ int main() {

cout << "Finished successfully!" << endl;
return 0;
}
}