From d2abdc21ef02a9f715aa7c03533671b3f968d24c Mon Sep 17 00:00:00 2001
From: Grigorii Drozdov <drozd013@umn.edu>
Date: Thu, 22 Apr 2021 22:35:16 -0500
Subject: [PATCH] fix: save exact particles state
---
apps/benchmarks/CNT/Config.cpp | 20 +++++++++++++-----
apps/benchmarks/CNT/InitializeCNTs.cpp | 29 ++++++++++++++++++++------
2 files changed, 38 insertions(+), 11 deletions(-)
diff --git a/apps/benchmarks/CNT/Config.cpp b/apps/benchmarks/CNT/Config.cpp
index 17e3b0896..aac837966 100644
--- a/apps/benchmarks/CNT/Config.cpp
+++ b/apps/benchmarks/CNT/Config.cpp
@@ -50,10 +50,11 @@ void saveConfig(const shared_ptr<data::ParticleStorage>& ps,
sb << static_cast<double> (p.getPosition()[1]);
sb << static_cast<double> (p.getPosition()[2]);
- Vec3 orient = p.getRotation().getMatrix() * Vec3(1,0,0);
- // Now orient has the shape (sin(th)cos(ph), sin(th)sin(ph), cos(th))
- sb << static_cast<double> (std::acos(orient[2]));
- sb << static_cast<double> (std::atan2(orient[1], orient[0]));
+ Quaternion q = p.getRotation().getQuaternion();
+ sb << static_cast<double> (q[0]);
+ sb << static_cast<double> (q[1]);
+ sb << static_cast<double> (q[2]);
+ sb << static_cast<double> (q[3]);
sb << static_cast<double> (p.getLinearVelocity()[0]);
sb << static_cast<double> (p.getLinearVelocity()[1]);
@@ -62,6 +63,15 @@ void saveConfig(const shared_ptr<data::ParticleStorage>& ps,
sb << static_cast<double> (p.getAngularVelocity()[0]);
sb << static_cast<double> (p.getAngularVelocity()[1]);
sb << static_cast<double> (p.getAngularVelocity()[2]);
+
+ sb << static_cast<double> (p.getOldForce()[0]);
+ sb << static_cast<double> (p.getOldForce()[1]);
+ sb << static_cast<double> (p.getOldForce()[2]);
+
+ sb << static_cast<double> (p.getOldTorque()[0]);
+ sb << static_cast<double> (p.getOldTorque()[1]);
+ sb << static_cast<double> (p.getOldTorque()[2]);
+
}
switch (io_mode)
@@ -74,7 +84,7 @@ void saveConfig(const shared_ptr<data::ParticleStorage>& ps,
WALBERLA_ROOT_SECTION()
{
uint_t dataSize = sizeof(mpi::RecvBuffer::ElementType) * rb.size();
- size_t size = rb.size() / 104;
+ size_t size = rb.size() / 168;
std::ofstream binfile;
binfile.open(filename + ".sav", std::ios::out | std::ios::binary);
binfile.write((char*) &size, sizeof(size_t));
diff --git a/apps/benchmarks/CNT/InitializeCNTs.cpp b/apps/benchmarks/CNT/InitializeCNTs.cpp
index ff9d29996..d37b3d460 100644
--- a/apps/benchmarks/CNT/InitializeCNTs.cpp
+++ b/apps/benchmarks/CNT/InitializeCNTs.cpp
@@ -164,27 +164,43 @@ int64_t loadCNTs(const std::string& filename,
double x;
double y;
double z;
- double theta;
- double phi;
+ double q0;
+ double q1;
+ double q2;
+ double q3;
double vx;
double vy;
double vz;
double wx;
double wy;
double wz;
+ double fx;
+ double fy;
+ double fz;
+ double tx;
+ double ty;
+ double tz;
binfile.read((char *) &sID, sizeof(int64_t));
binfile.read((char *) &cID, sizeof(int64_t));
binfile.read((char *) &x, sizeof(double));
binfile.read((char *) &y, sizeof(double));
binfile.read((char *) &z, sizeof(double));
- binfile.read((char *) &theta, sizeof(double));
- binfile.read((char *) &phi, sizeof(double));
+ binfile.read((char *) &q0, sizeof(double));
+ binfile.read((char *) &q1, sizeof(double));
+ binfile.read((char *) &q2, sizeof(double));
+ binfile.read((char *) &q3, sizeof(double));
binfile.read((char *) &vx, sizeof(double));
binfile.read((char *) &vy, sizeof(double));
binfile.read((char *) &vz, sizeof(double));
binfile.read((char *) &wx, sizeof(double));
binfile.read((char *) &wy, sizeof(double));
binfile.read((char *) &wz, sizeof(double));
+ binfile.read((char *) &fx, sizeof(double));
+ binfile.read((char *) &fy, sizeof(double));
+ binfile.read((char *) &fz, sizeof(double));
+ binfile.read((char *) &tx, sizeof(double));
+ binfile.read((char *) &ty, sizeof(double));
+ binfile.read((char *) &tz, sizeof(double));
Vec3 pos;
pos[0] = real_c(x);
pos[1] = real_c(y);
@@ -198,10 +214,11 @@ int64_t loadCNTs(const std::string& filename,
sp.setInteractionRadius(kernel::cnt::outer_radius);
sp.setSegmentID(sID);
sp.setClusterID(cID);
- sp.getRotationRef().rotate(Vec3(0_r, 1_r, 0_r), -0.5_r * math::pi + real_c(theta));
- sp.getRotationRef().rotate(Vec3(0_r, 0_r, 1_r), real_c(phi));
+ sp.setRotation(Rot3(Quaternion(q0, q1, q2, q3)));
sp.setLinearVelocity(Vec3(real_c(vx), real_c(vy), real_c(vz)));
sp.setAngularVelocity(Vec3(real_c(wx), real_c(wy), real_c(wz)));
+ sp.setOldForce(Vec3(real_c(fx), real_c(fy), real_c(fz)));
+ sp.setOldTorque(Vec3(real_c(tx), real_c(ty), real_c(tz)));
numParticles++;
}
}
--
GitLab