diff --git a/grains_in_c++/3D_grain_structure/datread.py b/grains_in_c++/3D_grain_structure/datread.py new file mode 100644 index 0000000..75e290b --- /dev/null +++ b/grains_in_c++/3D_grain_structure/datread.py @@ -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() diff --git a/grains_in_c++/3D_grain_structure/grain_maker.cc b/grains_in_c++/3D_grain_structure/grain_maker.cc index dd42348..7af7560 100644 --- a/grains_in_c++/3D_grain_structure/grain_maker.cc +++ b/grains_in_c++/3D_grain_structure/grain_maker.cc @@ -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 @@ -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; @@ -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) { @@ -69,25 +84,20 @@ int main() { } } } - + // --- KNOCK OUT GRAINS --- cout << "Knocking out invalid grains..." << endl; vector 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 @@ -106,12 +116,17 @@ int main() { // --- INITIALIZE ORDER PARAMETERS (PHI) --- vector 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; } } } @@ -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 Phi_new = Phi; vector 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) { @@ -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 --- @@ -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) { @@ -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 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 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()) { @@ -300,14 +319,12 @@ int main() { // 2. grains.dat (This must perfectly match the 4D interior space for your // other C++ script) - vector 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 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)]; } } } @@ -322,4 +339,4 @@ int main() { cout << "Finished successfully!" << endl; return 0; -} \ No newline at end of file +}