Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F72635488
lb_bonds_clusterwise.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Tue, Jul 16, 07:44
Size
23 KB
Mime Type
text/x-c
Expires
Thu, Jul 18, 07:44 (1 d, 21 h)
Engine
blob
Format
Raw Data
Handle
19081134
Attached To
rDLMA Diffusion limited mixed aggregation
lb_bonds_clusterwise.cpp
View Options
#include <post.h>
namespace post_p
{
void postprocessing::print_mem_usage()
{
getrusage(RUSAGE_SELF, &myusage);
std::cout<<" = "<<(1. * myusage.ru_maxrss)/1e3<<std::endl;
}
void postprocessing::dump_lb_bonds_for_cluster_via_invA(char *filename)
{
//std::cout<<"1 new"<<std::endl;
//getrusage(RUSAGE_SELF, &myusage);
//baseline = myusage.ru_maxrss;
//std::cout<<"on entering dump function = "<<(1. * baseline)/1e3<<std::endl;;
determine_LB_bonds_clusterwise(filename);
//std::cout<<"2"<<std::endl;
}
void postprocessing::calc_total_bonds()
{
total_num_bonds = 0;
for (int i = 0; i < numParticles(); i++)
total_num_bonds += numAttachments(i);
total_num_bonds /= 2;
}
bool postprocessing::check_if_particles_placed()
{
bool all_particles_placed = true;
for (int i = unf_checkpoint; i < numParticles(); i++) {
//all_particles_placed = all_particles_placed && is_placed(i);
if (is_placed(i) == false)
return false;
}
//std::cout<<"bool values = "<<all_particles_placed<<std::endl;
return true;
}
void postprocessing::distort_positions()
{
for (int i = 0; i < numParticles(); i++){
for (int axis = 0; axis < dim(); axis++){
//std::cout<<"old position = "<<pos(i,axis)<<"\t";
pos(i,axis) += ((0.1 * dis(generator))-0.05);
//std::cout<<"new position = "<<pos(i,axis)<<"\t";
}
//std::cout<<"\n";
}
//std::cout<<"distorted"<<std::endl;
}
void postprocessing::set_lbb_params(char *filename)
{
reset_unfolding_params();
build_bond_map();
reset_bond_map(true);
total_lbp = 0;
ref_pos.resize(dim());
totalClusters_ = 0;
f_lbp = fopen(filename, "w");
fprintf(f_lbp,"clusterNumber,bond,");
for (int axis = 0; axis < dim(); axis++)
fprintf(f_lbp,"bond_component_x%d,",axis);
fprintf(f_lbp,"bond_length\n");
/*if (!A.isCompressed()){
A.makeCompressed();
}*/
lbp_tolerance = 1e-6;
}
void postprocessing::init_lbb_unfolding()
{
num_bonds_for_cluster = 0;
num_particles_for_cluster = 0;
unique_bonds.clear();
index_to_particles.clear();
particles_to_index.clear();
for (int i = 0; i < numParticles(); i++){
if (is_placed(i) == false){
particles_to_index[i] = 0;
index_to_particles[0] = i;
num_particles_for_cluster += 1;
unfold_for_clusterwise(i,i);
break;
}
}
//std::cout<<"num_particles_for_cluster = "<<num_particles_for_cluster<<"\t num_bonds = "<<num_bonds_for_cluster<<std::endl;
}
void postprocessing::init_lbb_unfolding_without_recursion()
{
//unf_cp_1 = std::chrono::steady_clock::now();
num_bonds_for_cluster = 0;
num_particles_for_cluster = 0;
unique_bonds.clear();
index_to_particles.clear();
particles_to_index.clear();
to_build_list.clear();
//unf_cp_2 = std::chrono::steady_clock::now();
//unf_time_1 += std::chrono::duration_cast<std::chrono::nanoseconds>(unf_cp_2 - unf_cp_1).count();
//unf_cp_3 = std::chrono::steady_clock::now();
for (int i = unf_checkpoint; i < numParticles(); i++){
if (is_placed(i) == false){
particles_to_index[i] = num_particles_for_cluster;
index_to_particles[num_particles_for_cluster] = i;
to_build_list[num_particles_for_cluster] = i;
num_particles_for_cluster++;
unfold_for_clusterwise_without_recursion(i);
is_placed(i) = true;
for (int temp_index = 1; temp_index < num_particles_for_cluster; temp_index++){
unfold_for_clusterwise_without_recursion(to_build_list[temp_index]);
}
unf_checkpoint = i;
break;
}
}
//unf_cp_4 = std::chrono::steady_clock::now();
//unf_time_2 += std::chrono::duration_cast<std::chrono::nanoseconds>(unf_cp_4 - unf_cp_3).count();
//std::cout<<"num_particles_for_cluster = "<<num_particles_for_cluster<<"\t num_bonds = "<<num_bonds_for_cluster<<std::endl;
}
void postprocessing::unfold_for_clusterwise_without_recursion(int i)
{
for (int att = 0; att < numAttachments(i); att++)
{
index_j = attachment(i, att);
if (attachments_placed(index_j) == false)
{
if (is_placed(index_j) == false){
particles_to_index[index_j] = num_particles_for_cluster;
index_to_particles[num_particles_for_cluster] = index_j;
is_placed(index_j) = true;
//to_build_list.push_back(index_j);
to_build_list[num_particles_for_cluster] = index_j;
num_particles_for_cluster++;
//std::cout<<"index_i = "<<i<<"\t index_j = "<<index_j<<std::endl;
}
i_map = particles_to_index[i];
j_map = particles_to_index[index_j];
//unique_bonds.push_back({i_map,j_map});
unique_bonds[num_bonds_for_cluster] = {i_map,j_map};
num_bonds_for_cluster += 1;
//std::cout<<"i = "<<i<<"\t j = "<<index_j<<std::endl;
}
}
attachments_placed(i) = true;
}
void postprocessing::unfold_for_clusterwise(const int prev, const int next)
{
for (int axis = 0; axis < dim(); axis++)
posDiff(axis) = get_periodic_image(pos(next,axis) - pos(prev,axis), axis);
for (int axis = 0; axis < dim(); axis++){
if (prev == next)
unfolded_coords(next,axis) = pos(prev,axis) + posDiff(axis);
else
unfolded_coords(next,axis) = unfolded_coords(prev,axis) + posDiff(axis);
}
is_placed(next) = true;
for (int att = 0; att < numAttachments(next); att++)
{
temp_next = attachment(next, att);
if (is_placed(temp_next) == false){
particles_to_index[temp_next] = num_particles_for_cluster;
index_to_particles[num_particles_for_cluster] = temp_next;
num_particles_for_cluster++;
/*if (particles_to_index.count(next) != 1){
std::cout<<"unfold 1"<<std::endl;
}
if (particles_to_index.count(temp_next) != 1){
std::cout<<"unfold 2"<<std::endl;
}*/
index_i = particles_to_index[next];
index_j = particles_to_index[temp_next];
if (std::find(unique_bonds.begin(), unique_bonds.end(), std::make_pair(index_j,index_i)) == unique_bonds.end()){
unique_bonds.push_back({index_i,index_j});
num_bonds_for_cluster += 1;
}
unfoldedAttachments(next) = unfoldedAttachments(next) + 1;
unfoldedAttachments(temp_next) = unfoldedAttachments(temp_next) + 1;
unfold_for_clusterwise(next, temp_next);
}
else {
index_i = particles_to_index[next];
index_j = particles_to_index[temp_next];
if (std::find(unique_bonds.begin(), unique_bonds.end(), std::make_pair(index_j,index_i)) == unique_bonds.end()){
unique_bonds.push_back({index_i,index_j});
num_bonds_for_cluster += 1;
}
for (int axis = 0; axis < dim(); axis++) {
posDiff(axis) = unfolded_coords(temp_next, axis) - unfolded_coords(next, axis);
if ((posDiff(axis) > halfL(axis)) || (posDiff(axis) < -halfL(axis))){
cluster_percolation(totalClusters_,axis) = 1;
load_bearing_paths_[axis]++;
}
}
}
}
}
void postprocessing::init_lbb_cluster_matrices()
{
modified_folded_x.resize(num_particles_for_cluster,dim());
copy_positions_for_cluster();
nnz = 2 * num_bonds_for_cluster + num_particles_for_cluster;
A.reserve(nnz);
A.resize((num_particles_for_cluster-1), (num_particles_for_cluster-1));
b.resize((num_particles_for_cluster-1));
old_b.resize((num_particles_for_cluster-1));
x.resize((num_particles_for_cluster-1));
bond_lengths_direction_wise.resize(num_bonds_for_cluster, dim());
bond_lengths.resize(num_bonds_for_cluster);
A.setZero();
build_A_for_cluster();
//A.makeCompressed();
}
void postprocessing::print_lbb_info()
{
max_length = bond_lengths.maxCoeff(&max_row);
if (max_length > lbp_tolerance){
switch_off_bonds({index_to_particles[unique_bonds[max_row].first], index_to_particles[unique_bonds[max_row].second]});
modify_A_matrix();
fprintf(f_lbp, "%d,%d-%d,",totalClusters_,index_to_particles[unique_bonds[max_row].first],index_to_particles[unique_bonds[max_row].second]);
for (int print_axis = 0; print_axis < dim(); print_axis++)
fprintf(f_lbp, "%lf,", sqrt(bond_lengths_direction_wise(max_row,print_axis)));
fprintf(f_lbp, "%lf\n", max_length);
n_lbb++;
//std::cout<<"num_c = "<<num_particles_for_cluster<<"\t n_lbb = "<<n_lbb<<std::endl;
}
}
void postprocessing::print_minimized_coords(char *filename, int filenum)
{
std::string file_test(filename);
/*std::vector<std::string> results;
results = split_string_by_delimiter(file_test, '.');
std::string full_filename = results[0] + "_" + std::to_string(filenum) + "." + results[1];
for (int i = 0; i < results.size(); i++)
std::cout<<"results "<<i<<" = "<<results[i]<<std::endl;*/
file_test = file_test.substr(0, file_test.size()-4);
std::string full_filename = file_test + "_" + std::to_string(filenum) + ".csv";
std::ofstream f_coords;
f_coords.open(full_filename);
f_coords << "id,";
for (int axis = 0; axis < dim(); axis++){
if (axis == (dim() - 1)){
f_coords<<"x"<<axis<<"\n";
}
else{
f_coords<<"x"<<axis<<",";
}
}
int p_index;
for (int i = 0; i < num_particles_for_cluster; i++){
//p_index = particles_to_index[i];
p_index = index_to_particles[i];
//std::cout<<"i = "<<i<<"\t p_index = "<<std::endl;
f_coords<<p_index<<",";
for (int axis = 0; axis < dim(); axis++){
if (axis == (dim()-1)){
//f_coords<<modified_folded_x(p_index,axis)<<"\n";
f_coords<<modified_folded_x(i,axis)<<"\n";
}
else{
//f_coords<<modified_folded_x(p_index,axis)<<",";
f_coords<<modified_folded_x(i,axis)<<",";
}
}
}
f_coords.close();
}
void postprocessing::print_cluster_vals(char *filename)
{
std::string file_test(filename);
file_test = file_test.substr(0, file_test.size()-4);
std::string full_filename = file_test + "_clusters.csv";
std::ofstream f_coords;
f_coords.open(full_filename, std::ios_base::app);
for (int i = 0; i < num_particles_for_cluster; i++)
f_coords<<to_build_list[i]<<","<<totalClusters_<<"\n";
f_coords.close();
}
void postprocessing::determine_LB_bonds_clusterwise(char *filename)
{
std::chrono::steady_clock::time_point cp_1;
std::chrono::steady_clock::time_point cp_2;
std::chrono::steady_clock::time_point cp_3;
std::chrono::steady_clock::time_point cp_4;
std::chrono::steady_clock::time_point cp_5;
std::chrono::steady_clock::time_point cp_6;
std::chrono::steady_clock::time_point cp_7;
std::chrono::steady_clock::time_point cp_8;
std::string file_test(filename);
file_test = file_test.substr(0, file_test.size()-4);
std::string full_filename = file_test + "_clusters.csv";
std::ofstream f_coords;
f_coords.open(full_filename);
f_coords << "id,clusterNumber\n";
f_coords.close();
//std::cout<<"on entering determine function";print_mem_usage();
double unf_time = 0.;
double A_time = 0.;
double invA_time = 0.;
double while_loop_time = 0.;
int min_cluster_size = 3;
set_lbb_params(filename);
int num_A_constructions = 0;
int num_invAb_ops = 0;
//std::cout<<"before to_build_list";print_mem_usage();
to_build_list.resize(numParticles());
//std::cout<<"after to_build_list";print_mem_usage();
calc_total_bonds();
//std::cout<<"before unique bonds";print_mem_usage();
unique_bonds.resize(total_num_bonds);
//std::cout<<"after unique bonds";print_mem_usage();
//std::cout<<"after init bullshit";print_mem_usage();
totalClusters_ = 0;
while (!check_if_particles_placed()) {
//cp_7 = std::chrono::steady_clock::now();
//cp_1 = std::chrono::steady_clock::now();
//std::cout<<"before unfolding";print_mem_usage();
init_lbb_unfolding_without_recursion();
//std::cout<<"after unfolding";print_mem_usage();
//init_lbb_unfolding();
//cp_2 = std::chrono::steady_clock::now();
//unf_time += std::chrono::duration_cast<std::chrono::nanoseconds>(cp_2 - cp_1).count();
if (num_particles_for_cluster >= min_cluster_size) {
//std::cout<<"cluster size = "<<num_particles_for_cluster<<"\t num bonds = "<<num_bonds_for_cluster<<"\t time = "<<(unf_time * 1e-9)<<std::endl;
//cp_3 = std::chrono::steady_clock::now();
//std::cout<<"before init lbb matrices";print_mem_usage();
init_lbb_cluster_matrices();
//std::cout<<"after init lbb matrices";print_mem_usage();
//num_A_constructions++;
//cp_4 = std::chrono::steady_clock::now();
//A_time += std::chrono::duration_cast<std::chrono::nanoseconds>(cp_4 - cp_3).count();
//cp_5 = std::chrono::steady_clock::now();
do {
//std::cout<<"before modify coords";print_mem_usage();
modify_coords_for_cluster();
//std::cout<<"after modify coords";print_mem_usage();
if (!A.isCompressed()){
A.makeCompressed();
//std::cout<<"test"<<std::endl;
}
//std::cout<<"before analyze pattern";print_mem_usage();
solver.analyzePattern(A);
solver.factorize(A);
//std::cout<<"after analyze pattern";print_mem_usage();
//std::cout<<"lbb = "<<n_lbb<<std::endl;
for (int axis = 0; axis < dim(); axis++){
b.setZero();
build_b_for_cluster(axis);
x = solver.solve(b);
num_invAb_ops++;
modify_coords_after_minimization_for_cluster(axis);
}
//std::cout<<"after solving pattern";print_mem_usage();
calculate_bond_lengths_for_cluster();
//print_minimized_coords(filename, n_lbb);
print_lbb_info();
} while(max_length > lbp_tolerance);
//cp_6 = std::chrono::steady_clock::now();
//invA_time += std::chrono::duration_cast<std::chrono::nanoseconds>(cp_6 - cp_5).count();
}
print_cluster_vals(filename);
totalClusters_++;
//cp_8 = std::chrono::steady_clock::now();
//while_loop_time += std::chrono::duration_cast<std::chrono::nanoseconds>(cp_8 - cp_7).count();
//std::cout<<"total clusters = "<<totalClusters_<<std::endl;
}
//std::cout<<"while loop time = "<<(while_loop_time * 1e-9)<<std::endl;
fclose(f_lbp);
//std::cout<<"unfolding time = "<<unf_time * 1e-9<<std::endl;
//std::cout<<"A time = "<<A_time * 1e-9<<std::endl;
//std::cout<<"invA time = "<<invA_time * 1e-9<<std::endl;
//std::cout<<"n_lbb = "<<n_lbb<<std::endl;
//std::cout<<"num_A_constructions = "<<num_A_constructions<<std::endl;
//std::cout<<"num_invAb_ops = "<<num_invAb_ops<<std::endl;
//std::cout<<"unf 1 = "<<unf_time_1 * 1e-9<<std::endl;
//std::cout<<"unf 2 = "<<unf_time_2 * 1e-9<<std::endl;
//std::cout<<"A1 module = "<<time_A1<<std::endl;
//std::cout<<"A2 module = "<<time_A2<<std::endl;
}
void postprocessing::copy_positions_for_cluster()
{
int index;
for (int i = 0; i < num_particles_for_cluster; i++){
for (int axis = 0; axis < dim(); axis++){
index = index_to_particles[i];
modified_folded_x(i,axis) = pos(index,axis);
modified_folded_x(i,axis) -= L(axis) * round(modified_folded_x(i,axis)/L(axis));
}
}
}
void postprocessing::modify_coords_for_cluster()
{
int index;
for (int axis = 0; axis < dim(); axis++){
ref_pos(axis) = modified_folded_x(num_particles_for_cluster-1, axis);
}
for (int i = 0; i < num_particles_for_cluster; i++){
for (int axis = 0; axis < dim(); axis++){
modified_folded_x(i,axis) = modified_folded_x(i,axis) - ref_pos(axis);
modified_folded_x(i,axis) -= L(axis) * round(modified_folded_x(i,axis)/L(axis));
//std::cout<<pos(i,axis)<<" ";
}
//std::cout<<"\n";
}
}
void postprocessing::modify_coords_after_minimization_for_cluster(int axis)
{
for (int i = 0; i < num_particles_for_cluster-1; i++)
modified_folded_x(i,axis) = (L(axis) * x(i)) + ref_pos(axis);
modified_folded_x(num_particles_for_cluster-1,axis) = ref_pos(axis);
for (int i = 0; i < num_particles_for_cluster; i++)
modified_folded_x(i,axis) -= L(axis) * round(modified_folded_x(i,axis)/L(axis));
}
void postprocessing::put_particles_back_in_box(int axis)
{
for (int i = 0; i < num_particles_for_cluster; i++)
modified_folded_x(i,axis) -= L(axis) * round(modified_folded_x(i,axis)/L(axis));
}
void postprocessing::build_A_for_cluster()
{
int i,j;
//int index_i, index_j;
for (int index = 0; index < num_bonds_for_cluster; index++){
//cp_A1 = std::chrono::steady_clock::now();
i = unique_bonds[index].first;
j = unique_bonds[index].second;
//cp_A2 = std::chrono::steady_clock::now();
//time_A1 += std::chrono::duration_cast<std::chrono::nanoseconds>(cp_A2 - cp_A1).count();
//index_i = index_to_particles[i];
//index_j = index_to_particles[j];
/*if (bond_map_status.count({index_i,index_j}) == 0){
std::cout<<"problem with A"<<std::endl;
}*/
//if (bond_map_status[{index_i,index_j}] == 1){
cp_A3 = std::chrono::steady_clock::now();
if (i == (num_particles_for_cluster-1)){
tripleList.push_back(Eigen::Triplet<double>(j,j,-1));
}
else if (j == (num_particles_for_cluster-1)){
tripleList.push_back(Eigen::Triplet<double>(i,i,-1));
}
else{
tripleList.push_back(Eigen::Triplet<double>(i,j,1));
tripleList.push_back(Eigen::Triplet<double>(j,i,1));
tripleList.push_back(Eigen::Triplet<double>(i,i,-1));
tripleList.push_back(Eigen::Triplet<double>(j,j,-1));
}
//cp_A4 = std::chrono::steady_clock::now();
//time_A2 += std::chrono::duration_cast<std::chrono::nanoseconds>(cp_A4 - cp_A3).count();
//}
}
A.setFromTriplets(tripleList.begin(), tripleList.end());
tripleList.clear();
}
void postprocessing::build_b_for_cluster(int axis)
{
int i,j;
int index_i, index_j;
double dx;
for (int index = 0; index < num_bonds_for_cluster; index++){
i = unique_bonds[index].first;
j = unique_bonds[index].second;
index_i = index_to_particles[i];
index_j = index_to_particles[j];
/*if (bond_map_status.count({index_i,index_j}) == 0){
std::cout<<"problem with b"<<std::endl;
}*/
if (bond_map_status[{index_i,index_j}] == 1){
dx = round((modified_folded_x(i,axis) - modified_folded_x(j,axis))/L(axis));
if (i == (num_particles_for_cluster-1)){
b(j) += dx;
}
else if (j == (num_particles_for_cluster-1)){
b(i) -= dx;
}
else {
b(i) -= dx;
b(j) += dx;
}
}
}
}
void postprocessing::modify_A_matrix()
{
a_i = unique_bonds[max_row].first;
a_j = unique_bonds[max_row].second;
if (a_i == (num_particles_for_cluster-1)){
A.coeffRef(a_j,a_j) += 1;
}
else if (a_j == (num_particles_for_cluster-1)){
A.coeffRef(a_i, a_i) += 1;
}
else{
A.coeffRef(a_i,a_j) = 0;
A.coeffRef(a_j,a_i) = 0;
A.coeffRef(a_i,a_i) += 1;
A.coeffRef(a_j,a_j) += 1;
}
}
void postprocessing::copy_b()
{
for (int i = 0; i < b.size(); i++)
old_b(i) = b(i);
}
void postprocessing::calculate_b_diff()
{
b_diff = b - old_b;
if (b_diff.norm() > 1e-4)
std::cout<<"b vector changes"<<std::endl;
}
void postprocessing::increase_stiffness()
{
for (int i = 0; i < b.size(); i++){
b(i) = 1. * b(i);
}
}
void postprocessing::calculate_bond_lengths_direction_wise_for_cluster(int axis)
{
//bond_distance.setZero();
int i,j;
int index_i, index_j;
for (int index = 0; index < num_bonds_for_cluster; index++){
i = unique_bonds[index].first;
j = unique_bonds[index].second;
index_i = index_to_particles[i];
index_j = index_to_particles[j];
if (bond_map_status[{index_i,index_j}] == 1){
bond_lengths_direction_wise(index, axis) = (modified_folded_x(i, axis) - modified_folded_x(j, axis));
bond_lengths_direction_wise(index, axis) -= (L(axis) * round(bond_lengths_direction_wise(index, axis)/L(axis)));
bond_lengths_direction_wise(index, axis) *= bond_lengths_direction_wise(index, axis);
}
else {
bond_lengths_direction_wise(index, axis) = 0.;
}
}
}
void postprocessing::calculate_bond_lengths_for_cluster()
{
bond_lengths_direction_wise.setZero();
bond_lengths.setZero();
for (int axis = 0; axis < dim(); axis++){
calculate_bond_lengths_direction_wise_for_cluster(axis);
}
for (int i = 0; i < num_bonds_for_cluster; i++){
for (int axis = 0; axis < dim(); axis++) {
bond_lengths(i) += bond_lengths_direction_wise(i,axis);
}
bond_lengths(i) = sqrt(bond_lengths(i));
}
}
void postprocessing::calculate_bond_lengths_for_cluster(int axis)
{
//bond_distance.setZero();
int i,j;
int index_i, index_j;
for (int index = 0; index < num_bonds_for_cluster; index++){
i = unique_bonds[index].first;
j = unique_bonds[index].second;
index_i = index_to_particles[i];
index_j = index_to_particles[j];
if (bond_map_status[{index_i,index_j}] == 1){
bond_lengths(index) = (modified_folded_x(i, axis) - modified_folded_x(j, axis));
bond_lengths(index) -= (L(axis) * round(bond_lengths(index)/L(axis)));
bond_lengths(index) = abs(bond_lengths(index));
}
else {
bond_lengths(index) = 0.;
}
//std::cout<<index_i<<"-"<<index_j<<" = "<<bond_lengths[index]<<std::endl;
}
//std::cout<<"-----------------------------------------------------"<<std::endl;
}
void postprocessing::dump_xyz_stepwise(char *filename, int file_num)
{
std::string ext_filename;
ext_filename = filename;
ext_filename += "_" + std::to_string(file_num) + ".csv";
std::ofstream test_f;
test_f.open(ext_filename);
for (int axis = 0; axis < dim(); axis++){
if (axis != (dim()-1))
test_f <<"x"+std::to_string(axis)+",";
else
test_f <<"x"+std::to_string(axis)+"\n";
}
for (int i = 0; i < numParticles(); i++){
for (int axis = 0; axis < dim(); axis++){
if (axis != (dim()-1))
test_f <<std::to_string(modified_folded_x(i,axis))+",";
else
test_f <<std::to_string(modified_folded_x(i,axis))+"\n";
}
}
test_f.close();
}
}
Event Timeline
Log In to Comment