commit 0e97bc9c31d39275196935a69d44fcd9029cbb35
parent a44793fa3b64c460c2cb675919424df4e888932b
Author: Ursula Rasthofer <urasthofer@ethz.ch>
Date: Mon, 6 Feb 2017 18:55:54 +0100
added RP-cluster
Diffstat:
5 files changed, 244 insertions(+), 2 deletions(-)
diff --git a/Makefile b/Makefile
@@ -7,6 +7,7 @@ sse ?= 0
align ?= 16
sundials ?= 0
eigen ?= 0
+pos ?= 0
sos = sg
LIBS =
@@ -51,6 +52,10 @@ ifeq "$(eigen)" "1"
CXXFLAGS += -I/usr/local/opt/eigen/include/eigen3
endif
+ifeq "$(pos)" "1"
+ CXXFLAGS += -D_POS_STATE_
+endif
+
VPATH := src
SRC = src/bubbleDynamics.cpp
HDR = $(wildcard *.h)
diff --git a/src/bubbleDynamics.cpp b/src/bubbleDynamics.cpp
@@ -43,6 +43,7 @@ KernelBase<Tinput,Trhs>* kernelFactory(ArgumentParser& parser, Tinput& U, Bubble
else
yeaBinary = GnuplotDump::ASCII;
+#ifndef _POS_STATE_
if (parser("-dynamics").asString("rp") == "rp")
{
RayleighPlesset<Tinput,Trhs>::IC(U, simConfig);
@@ -64,9 +65,28 @@ KernelBase<Tinput,Trhs>* kernelFactory(ArgumentParser& parser, Tinput& U, Bubble
KMCluster_TY<Tinput,Trhs>::IC(U, simConfig);
return new KMCluster_TY<Tinput,Trhs>(U.size(), yeaBinary);
}
+ else if (parser("-dynamics").asString("rp") == "kmcluster-fc")
+ {
+ KMCluster_FC<Tinput,Trhs>::IC(U, simConfig);
+ return new KMCluster_FC<Tinput,Trhs>(U.size(), yeaBinary);
+ }
+ else if (parser("-dynamics").asString("rp") == "rpcluster")
+ {
+ RPCluster<Tinput,Trhs>::IC(U, simConfig);
+ return new RPCluster<Tinput,Trhs>(U.size(), yeaBinary);
+ }
#endif
else
return nullptr;
+#else
+ if (parser("-dynamics").asString("rp") == "kmcluster-positions")
+ {
+ KMClusterPositions_D<Tinput,Trhs>::IC(U, simConfig);
+ return new KMClusterPositions_D<Tinput,Trhs>(U.size(), yeaBinary);
+ }
+ else
+ return nullptr;
+#endif
}
template <typename Tinput, typename Trhs=Tinput>
@@ -92,7 +112,11 @@ int main(int argc, const char** argv)
BubbleData simConfig(parser);
// set up solution vector
+#if (_POS_STATE_ && _USE_EIGEN_)
+ using vec_t = StateVector<State8>;
+#else
using vec_t = StateVector<State2>;
+#endif
vec_t U(simConfig.Nbubbles);
// Select stepper and kernel
diff --git a/src/kernels.h b/src/kernels.h
@@ -10,7 +10,12 @@
#include "kernels/RayleighPlesset_HBGL.h"
#include "kernels/KellerMiksis.h"
#ifdef _USE_EIGEN_
+#include "kernels/RPCluster.h"
#include "kernels/KMCluster_TY.h"
+#include "kernels/KMCluster_FC.h"
+#ifdef _POS_STATE_
+#include "kernels/KMClusterPositions_D.h"
+#endif
#endif
#endif /* KERNELS_H_KY4W86QV */
diff --git a/src/kernels/KMCluster_TY.h b/src/kernels/KMCluster_TY.h
@@ -12,8 +12,10 @@
#include <cmath>
#include <vector>
-#include <Eigen/Core>
-#include <Eigen/Dense>
+#include "/cluster/home/ursular/research/codes/eigen/Eigen/Core"
+#include "/cluster/home/ursular/research/codes/eigen/Eigen/Dense"
+//#include <Eigen/Core>
+//#include <Eigen/Dense>
#include "KernelBase.h"
#include "BubbleBase.h"
diff --git a/src/kernels/RPCluster.h b/src/kernels/RPCluster.h
@@ -0,0 +1,206 @@
+/* File: RPCluster.h */
+/* Date: Sun Feb 05 22:10:34 2017 */
+/* Author: Ursula Rasthofer */
+/* Tag: Bubble cluster based on Rayleigh-Plesset */
+/* Copyright 2017 ETH Zurich. All Rights Reserved. */
+#ifndef RPCLUSTER_H_67YUOP08
+#define RPCLUSTER_H_67YUOP08
+
+#include <cassert>
+#include <cstdlib>
+#include <fstream>
+#include <cmath>
+#include <vector>
+
+#include "/cluster/home/ursular/research/codes/eigen/Eigen/Core"
+#include "/cluster/home/ursular/research/codes/eigen/Eigen/Dense"
+//#include <Eigen/Core>
+//#include <Eigen/Dense>
+
+#include "KernelBase.h"
+#include "BubbleBase.h"
+#include "GnuplotDump.h"
+
+
+template <typename Tinput, typename Trhs=Tinput>
+class RPCluster : public KernelBase<Tinput,Trhs>
+{
+ const size_t _N;
+ GnuplotDump* m_dumper;
+
+public:
+ RPCluster(const size_t N=0, const int format=GnuplotDump::BIN) : _N(N) { m_dumper = new GnuplotDump("out", format); }
+ virtual ~RPCluster() { delete m_dumper;}
+
+ static void IC(Tinput& U, BubbleData& data)
+ {
+ assert(Tinput::DataType::SIZE == 2);
+ assert(data.Nbubbles == U.size());
+ const size_t N = data.Nbubbles;
+
+ // read IC from file
+ std::ifstream infile(data.filename);
+ size_t k;
+ infile >> k;
+ if (N != k)
+ {
+ std::cerr << "RPCluster: Allocated " << N << " bubbles."<< std::endl;
+ std::cerr << " File " << data.filename << " expects " << k << " bubbles."<< std::endl;
+ abort();
+ }
+
+ for (size_t i = 0; i < N; ++i)
+ {
+ if (infile.good())
+ {
+ Bubble& b = data.coords[i];
+ infile >> b.pos[0] >> b.pos[1] >> b.pos[2];
+ infile >> data.R0[i];
+ infile >> data.Rdot0[i];
+
+ State2& IC = U[i];
+ IC[0] = data.R0[i];
+ IC[1] = data.Rdot0[i];
+ }
+ else
+ {
+ std::cerr << "RPCluster: Input file has incorrect format" << std::endl;
+ abort();
+ }
+ }
+
+ for (int i = 0; i < N; ++i)
+ {
+ const Bubble& bi = data.coords[i];
+ for (int j = 0; j < N; ++j)
+ {
+ const Bubble& bj = data.coords[j];
+ data.Dij[j + i*N] = bi.distance(bj);
+ }
+ }
+
+#ifndef NDEBUG
+ for (size_t i = 0; i < N; ++i)
+ assert(0 == data.Dij[i + i*N]);
+#endif /* NDEBUG */
+ }
+
+ virtual void compute(const Tinput& U, Trhs& rhs, Real const t, void const* const data=nullptr)
+ {
+ const BubbleData& bd = *(BubbleData const* const)data;
+
+ assert(_N == U.size());
+ Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic> A(_N,_N);
+ Eigen::Matrix<Real, Eigen::Dynamic, 1> b(_N);
+
+ for (size_t i = 0; i < _N; ++i)
+ {
+ const State2& Ui = U[i];
+ // fill row i of A
+ Real bnbr = 0;
+ for (size_t j = 0; j < _N; ++j)
+ {
+ const State2& Uj = U[j];
+ if (i == j)
+ A(i,j) = Ui[0];
+ else
+ {
+ const Real UDij = Uj[0]/bd.Dij[j + i*_N];
+ A(i,j) = Uj[0]*UDij;
+ bnbr += Uj[1]*Uj[1]*UDij;
+ }
+ }
+
+ // compute bi
+ const Real pstat = bd.p0 - bd.pv;
+ const Real pbs = (pstat + static_cast<Real>(2.0)*bd.sigma/bd.R0[i])*std::pow(bd.R0[i]/Ui[0], 3.0*bd.gamma);
+
+ b(i) = 0.0;
+ b(i) -= static_cast<Real>(1.5)*Ui[1]*Ui[1];
+ b(i) -= bd.rInv*(pstat + bd.pext(t));
+ b(i) += bd.rInv*(pbs - static_cast<Real>(2.0)*bd.sigma/Ui[0]) - static_cast<Real>(4.0)*bd.nuL*Ui[1]/Ui[0];
+ b(i) -= static_cast<Real>(2.0)*bnbr;
+
+ assert(!isnan(b(i)));
+ }
+
+ /* std::cout << A << std::endl << std::endl; */
+ //Eigen::Matrix<Real, Eigen::Dynamic, 1> Rddot = A.fullPivHouseholderQr().solve(b);
+ // this is notably faster and sufficiently accurate (rasthofer Feb 3, 2017)
+ Eigen::Matrix<Real, Eigen::Dynamic, 1> Rddot = A.colPivHouseholderQr().solve(b);
+
+ for (size_t i = 0; i < _N; ++i)
+ {
+ const State2& Ui = U[i];
+ State2& ri = rhs[i];
+ ri[0] = Ui[1];
+ ri[1] = Rddot(i);
+ }
+ }
+
+ virtual void write(const size_t step, const Real t, const Real dt, const Tinput& U, const void * const data=nullptr)
+ {
+ static size_t callCount = 0;
+
+ const BubbleData& bd = *(BubbleData const* const)data;
+
+ // double pb[3];
+ // const int id[3] = {0, 7, 1};
+ // for (int i = 0; i < 3; ++i)
+ // pb[i] = bd.pBubble(U, id[i]);
+
+ // // class 1, 2, 3, 4, 5, 6
+ // // const Real odat[6] = {U[0][0], U[1][0], U[4][0], U[2][0], U[5][0], U[3][0]};
+ // const Real odat[7] = {U[0][0], U[7][0], U[1][0], pb[0], pb[1], pb[2], bd.pext(t)};
+
+ const size_t N = 2*_N+2;
+ double out[N];
+ double Vg = 0.0;
+ for (int i = 0; i < _N; ++i)
+ {
+ const double r = U[i][0];
+ out[i] = r;
+ out[i+_N] = bd.pBubble(U, i);
+ Vg += 4.0/3.0*r*r*r*M_PI;
+ }
+ out[2*_N] = bd.pext(t);
+ out[2*_N+1] = Vg;
+
+ m_dumper->write(step, t, dt, &out[0], N);
+
+ if (bd.bVTK)
+ {
+ // dump bubble vtk
+ std::stringstream streamer;
+ streamer<< "bubbles-vtk_" << setfill('0') << setw(8) << callCount++ << ".vtk";
+ std::ofstream vtk(streamer.str().c_str());
+ vtk.setf(std::ios::scientific, std::ios::floatfield);
+
+ vtk << "# vtk DataFile Version 2.0\n";
+ vtk << "Unstructured grid legacy vtk file with point scalar data\n";
+ vtk << "ASCII\n" << std::endl;
+
+ vtk << "DATASET UNSTRUCTURED_GRID\n";
+ vtk << "POINTS " << bd.Nbubbles << " double" << std::endl;
+ for (int i = 0; i < bd.Nbubbles; ++i)
+ vtk << bd.coords[i].pos[0] << " " << bd.coords[i].pos[1] << " " << bd.coords[i].pos[2] << std::endl;
+ vtk << std::endl;
+
+ vtk << "POINT_DATA " << bd.Nbubbles << std::endl;
+ vtk << "SCALARS radius double\n";
+ vtk << "LOOKUP_TABLE default" << std::endl;
+ for (int i = 0; i < bd.Nbubbles; ++i)
+ vtk << U[i][0] << std::endl;
+ vtk << std::endl;
+
+ // vtk << "POINT_DATA " << bd.Nbubbles << std::endl;
+ vtk << "SCALARS velocity double\n";
+ vtk << "LOOKUP_TABLE default" << std::endl;
+ for (int i = 0; i < bd.Nbubbles; ++i)
+ vtk << U[i][1] << std::endl;
+ vtk << std::endl;
+ }
+ }
+};
+
+#endif /* RPCLUSTER_H_67YUOP08 */
