commit 7aa18619aa9c531bfc74550d1dbe30c34c2b22e4
parent 0e97bc9c31d39275196935a69d44fcd9029cbb35
Author: Ursula Rasthofer <urasthofer@ethz.ch>
Date: Tue, 7 Feb 2017 11:32:11 +0100
bubble cluster with particle translation according to Doinikov 2004
KM and RP version
some minor clean-up of the code will follow
Diffstat:
4 files changed, 701 insertions(+), 1 deletion(-)
diff --git a/src/bubbleDynamics.cpp b/src/bubbleDynamics.cpp
@@ -79,7 +79,12 @@ KernelBase<Tinput,Trhs>* kernelFactory(ArgumentParser& parser, Tinput& U, Bubble
else
return nullptr;
#else
- if (parser("-dynamics").asString("rp") == "kmcluster-positions")
+ if (parser("-dynamics").asString("rp") == "rpcluster-positions")
+ {
+ RPClusterPositions_D<Tinput,Trhs>::IC(U, simConfig);
+ return new RPClusterPositions_D<Tinput,Trhs>(U.size(), yeaBinary);
+ }
+ else if (parser("-dynamics").asString("rp") == "kmcluster-positions")
{
KMClusterPositions_D<Tinput,Trhs>::IC(U, simConfig);
return new KMClusterPositions_D<Tinput,Trhs>(U.size(), yeaBinary);
diff --git a/src/kernels.h b/src/kernels.h
@@ -14,6 +14,7 @@
#include "kernels/KMCluster_TY.h"
#include "kernels/KMCluster_FC.h"
#ifdef _POS_STATE_
+#include "kernels/RPClusterPositions_D.h"
#include "kernels/KMClusterPositions_D.h"
#endif
#endif
diff --git a/src/kernels/KMClusterPositions_D.h b/src/kernels/KMClusterPositions_D.h
@@ -0,0 +1,350 @@
+/* File: KMClusterPositions_D.h */
+/* Date: Sun Feb 5 13:24:54 2017 */
+/* Author: Ursula Rasthofer */
+/* Tag: Bubble cluster based on Keller-Miksis */
+/* including translation of bubbles following Doinikov (2004) */
+/* Copyright 2017 ETH Zurich. All Rights Reserved. */
+#ifndef BUBBLECLUSTERPOS_D_H_GHBNU569
+#define BUBBLECLUSTERPOS_D_H_GHBNU569
+
+#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 KMClusterPositions_D : public KernelBase<Tinput,Trhs>
+{
+ const size_t _N;
+ GnuplotDump* m_dumper;
+
+public:
+ KMClusterPositions_D(const size_t N=0, const int format=GnuplotDump::BIN) : _N(N) { m_dumper = new GnuplotDump("out", format); }
+ virtual ~KMClusterPositions_D() { delete m_dumper;}
+
+ static void IC(Tinput& U, BubbleData& data)
+ {
+ assert(Tinput::DataType::SIZE == 8);
+ 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 << "KMClusterPositions_D: 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];
+
+ State8& IC = U[i];
+ IC[0] = data.R0[i];
+ IC[1] = data.Rdot0[i];
+ IC[2] = b.pos[0];
+ IC[3] = 0.0;
+ IC[4] = b.pos[1];
+ IC[5] = 0.0;
+ IC[6] = b.pos[2];
+ IC[7] = 0.0;
+ }
+ else
+ {
+ std::cerr << "KMClusterPositions_D: 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());
+ const int size = 4*_N;
+ Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic> A(size,size);
+ Eigen::Matrix<Real, Eigen::Dynamic, 1> b(size);
+
+ const Real clInv = 1.0/bd.cl;
+
+ for (size_t i = 0; i < _N; ++i)
+ {
+ const State8& Ui = U[i];
+ const Bubble bi = {Ui[2], Ui[4], Ui[6]};
+ // fill row i of A
+ Real bnbr = 0;
+ std::vector<Real> bnbrp(3,0.0);
+ for (size_t j = 0; j < _N; ++j)
+ {
+ const State8& Uj = U[j];
+ if (i == j)
+ {
+ A(i,j) = (static_cast<Real>(1) - Ui[1]*clInv)*Ui[0] + static_cast<Real>(4.0)*bd.nuL*clInv;
+ for (std::size_t rr=0; rr<3; rr++)
+ {
+ const int ip = _N+3*i+rr;
+ A(i,ip) = static_cast<Real>(0.0);
+ A(ip,ip) = Ui[0]/static_cast<Real>(3.0);
+ A(ip,_N+3*i+(rr+1)%3) = static_cast<Real>(0.0);
+ A(ip,_N+3*i+(rr+2)%3) = static_cast<Real>(0.0);
+ A(ip,j) = static_cast<Real>(0.0);
+ }
+ }
+ else
+ {
+ // compute distance between bubbles
+ const Bubble bj = {Uj[2], Uj[4], Uj[6]};
+ const Real Dij = bi.distance(bj);
+
+ if (Dij < (Ui[0]+Uj[0]))
+ {
+ std::cerr << "Colliding bubbles " << Dij << " " << Ui[0] << " " << Uj[0] << std::endl;
+ abort();
+ }
+
+ const Real UDij = Uj[0]/Dij;
+ // RR-block
+ A(i,j) = Uj[0]*UDij;
+ bnbr += Uj[1]*Uj[1]*UDij; // keep in mind: multiplication by factor 2 below!!
+ // Rp/pp/pR-block
+ const Real Uj2Dij3 = UDij*UDij/Dij;
+ // factors for matrix
+ //
+ // Rj^3
+ // ---------
+ // 2 dij^3
+ const Real RpUDij3 = static_cast<Real>(0.5)*Uj2Dij3*Uj[0];
+ //
+ // - Ri Rj^2
+ // ----------
+ // dij^3
+ const Real pRUDij3 = static_cast<Real>(-1.0)*Uj2Dij3*Ui[0];
+ //
+ // Ri Rj^3
+ // ---------
+ // 2 dij^3
+ const Real ppUDij3 = RpUDij3*Ui[0];
+ //
+ // - 3 Ri Rj^3
+ // ------------
+ // 2 dij^5
+ const Real ppUDij5 = static_cast<Real>(-3.0)*ppUDij3/(Dij*Dij);
+ // factors for rhs
+ //
+ // radius part
+ //
+ // Rj^2 Rdj
+ // --------- (comment 1/4 due to multiplication by factor 2 below)
+ // 2 dij^3
+ const Real rhsR1 = static_cast<Real>(0.25)*Uj2Dij3*Uj[1];
+ //
+ // Rj^3
+ // --------- (comment 1/8 due to multiplication by factor 2 below)
+ // 4 dij^3
+ const Real rhsR2 = static_cast<Real>(0.125)*Uj2Dij3*Uj[0];
+ //
+ // 3 Rj^3
+ // --------- (comment 3/8 due to multiplication by factor 2 below)
+ // 4 dij^5
+ const Real velpro = (Uj[2] - Ui[2])*Uj[3] + (Uj[4] - Ui[4])*Uj[5] + (Uj[6] - Ui[6])*Uj[7];
+ const Real rhsR3 = static_cast<Real>(3.0)*rhsR2*velpro/(Dij*Dij);
+ //
+ // position part
+ //
+ // Rj Rdj
+ // -------- (2 Ri Rdj + Rdi Rj)
+ // dij^3
+ const Real rhsp1 = (static_cast<Real>(2.0)*Ui[0]*Uj[1]+Ui[1]*Uj[0])*UDij*Uj[1]/(Dij*Dij);
+ //
+ // Rj^2
+ // -------- (Rdi Rj + 5 Ri Rdj)
+ // 2 dij^3
+ const Real rhsp2 = static_cast<Real>(0.5)*Uj2Dij3*(Ui[1]*Uj[0]+static_cast<Real>(5.0)*Ui[0]*Uj[1]);
+ //
+ // 3 Rj^2
+ // -------- (Rdi Rj + 5 Ri Rdj)
+ // 2 dij^5
+ const Real rhsp3 = static_cast<Real>(3.0)*rhsp2/(Dij*Dij) * velpro; // consider opposite sign of velpro below (substraction instead of addition)
+ for (std::size_t rr=0; rr<3; rr++)
+ {
+ const int jp = _N+3*j+rr;
+ const int ip = _N+3*i+rr;
+ const Real pos_ij = Ui[rr*2+2] - Uj[rr*2+2];
+ // Rp-block
+ A(i,jp) = RpUDij3*pos_ij;
+ // pR-block
+ A(ip,j) = pRUDij3*pos_ij;
+ // pp-block
+ A(ip,jp) = ppUDij3;
+ A(ip,jp) += ppUDij5*pos_ij*pos_ij;
+ const int ia = (rr+1)%3;
+ const int ib = (rr+2)%3;
+ A(ip,_N+3*j+ia) = ppUDij5*pos_ij*(Ui[ia*2+2] - Uj[ia*2+2]);
+ A(ip,_N+3*j+ib) = ppUDij5*pos_ij*(Ui[ib*2+2] - Uj[ib*2+2]);
+
+ // rhs contributions
+ Real vel_ij = Ui[rr*2+3] + static_cast<Real>(5.0) * Uj[rr*2+3];
+ bnbr += rhsR1*pos_ij*vel_ij;
+ vel_ij = Ui[rr*2+3] + static_cast<Real>(2.0) * Uj[rr*2+3];
+ bnbr -= rhsR2*Uj[rr*2+3]*vel_ij;
+ bnbr -= rhsR3*pos_ij*vel_ij;
+ bnbrp[rr] += rhsp1*pos_ij;
+ bnbrp[rr] -= rhsp2*Uj[rr*2+3];
+ bnbrp[rr] -= rhsp3*pos_ij; //TODO: save multiplication here precomputing (rhsp1-rhsp3)
+ }
+ }
+ }
+
+ // compute bi
+ const Real pstat = bd.p0 - bd.pv;
+ const Real Rdc = Ui[1]*clInv;
+ 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) = (static_cast<Real>(0.5)*Rdc - static_cast<Real>(1.5))*Ui[1]*Ui[1];
+ b(i) -= bd.rInv*(static_cast<Real>(1.0) + Rdc)*(pstat + bd.pext(t + Ui[0]*clInv));
+ // b(i) -= bd.rInv*(static_cast<Real>(1.0) + Rdc)*(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) += bd.rInv*static_cast<Real>(1.0 - 3.0*bd.gamma)*Rdc*pbs;
+ b(i) -= bd.rInv*clInv*Ui[0]*bd.dpext(t + Ui[0]*clInv);
+ // b(i) -= bd.rInv*Rdc*bd.dpext(t);
+ // coupling terms
+ b(i) -= static_cast<Real>(2.0)*bnbr;
+ // translation term
+ b(i) += static_cast<Real>(0.25)*(Ui[3]*Ui[3] + Ui[5]*Ui[5] + Ui[7]*Ui[7]);
+
+ // pressure part
+ for (std::size_t rr=0; rr<3; rr++)
+ {
+ const int ip = _N+3*i+rr;
+ // coupling terms
+ b(ip) = bnbrp[rr];
+ // std rhs
+ b(ip) -= Ui[1]*Ui[rr*2+3];
+ // primary Bjerknes force and viscous drag neglected
+ }
+ assert(!isnan(b(i)));
+ }
+
+ //std::cout << A << std::endl << std::endl;
+ //std::cout << b << 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 State8& Ui = U[i];
+ State8& ri = rhs[i];
+ ri[0] = Ui[1];
+ ri[1] = Rddot(i);
+ ri[2] = Ui[3];
+ ri[3] = Rddot(_N+3*i);
+ ri[4] = Ui[5];
+ ri[5] = Rddot(_N+3*i+1);
+ ri[6] = Ui[7];
+ ri[7] = Rddot(_N+3*i+2);
+ }
+ }
+
+ 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;
+
+ const size_t N = 5*_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);
+ out[i+2*_N] = U[i][2];
+ out[i+3*_N] = U[i][4];
+ out[i+4*_N] = U[i][6];
+ Vg += 4.0/3.0*r*r*r*M_PI;
+ }
+ out[5*_N] = bd.pext(t);
+ out[5*_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 << U[i][2] << " " << U[i][4] << " " << U[i][6] << 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;
+
+ vtk << "VECTORS transvel double\n";
+ vtk << "LOOKUP_TABLE default" << std::endl;
+ for (int i = 0; i < bd.Nbubbles; ++i)
+ vtk << U[i][3] << " " << U[i][5] << " " << U[i][7] << std::endl;
+ vtk << std::endl;
+ }
+ }
+};
+
+#endif /* BUBBLECLUSTERPOS_D_H_GHBNU569 */
diff --git a/src/kernels/RPClusterPositions_D.h b/src/kernels/RPClusterPositions_D.h
@@ -0,0 +1,344 @@
+/* File: RPClusterPositions_D.h */
+/* Date: Tue Feb 7 09:26:31 2017 */
+/* Author: Ursula Rasthofer */
+/* Tag: Bubble cluster based on Rayleigh-Plesset */
+/* including translation of bubbles following Doinikov (2004) */
+/* Copyright 2017 ETH Zurich. All Rights Reserved. */
+#ifndef RPCLUSTERPOS_D_H_UIJM88YU
+#define RPCLUSTERPOS_D_H_UIJM88YU
+
+#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 RPClusterPositions_D : public KernelBase<Tinput,Trhs>
+{
+ const size_t _N;
+ GnuplotDump* m_dumper;
+
+public:
+ RPClusterPositions_D(const size_t N=0, const int format=GnuplotDump::BIN) : _N(N) { m_dumper = new GnuplotDump("out", format); }
+ virtual ~RPClusterPositions_D() { delete m_dumper;}
+
+ static void IC(Tinput& U, BubbleData& data)
+ {
+ assert(Tinput::DataType::SIZE == 8);
+ 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 << "RPClusterPositions_D: 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];
+
+ State8& IC = U[i];
+ IC[0] = data.R0[i];
+ IC[1] = data.Rdot0[i];
+ IC[2] = b.pos[0];
+ IC[3] = 0.0;
+ IC[4] = b.pos[1];
+ IC[5] = 0.0;
+ IC[6] = b.pos[2];
+ IC[7] = 0.0;
+ }
+ else
+ {
+ std::cerr << "RPClusterPositions_D: 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());
+ const int size = 4*_N;
+ Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic> A(size,size);
+ Eigen::Matrix<Real, Eigen::Dynamic, 1> b(size);
+
+ for (size_t i = 0; i < _N; ++i)
+ {
+ const State8& Ui = U[i];
+ const Bubble bi = {Ui[2], Ui[4], Ui[6]};
+ // fill row i of A
+ Real bnbr = 0;
+ std::vector<Real> bnbrp(3,0.0);
+ for (size_t j = 0; j < _N; ++j)
+ {
+ const State8& Uj = U[j];
+ if (i == j)
+ {
+ A(i,j) = Ui[0];
+ for (std::size_t rr=0; rr<3; rr++)
+ {
+ const int ip = _N+3*i+rr;
+ A(i,ip) = static_cast<Real>(0.0);
+ A(ip,ip) = Ui[0]/static_cast<Real>(3.0);
+ A(ip,_N+3*i+(rr+1)%3) = static_cast<Real>(0.0);
+ A(ip,_N+3*i+(rr+2)%3) = static_cast<Real>(0.0);
+ A(ip,j) = static_cast<Real>(0.0);
+ }
+ }
+ else
+ {
+ // compute distance between bubbles
+ const Bubble bj = {Uj[2], Uj[4], Uj[6]};
+ const Real Dij = bi.distance(bj);
+
+ if (Dij < (Ui[0]+Uj[0]))
+ {
+ std::cerr << "Colliding bubbles " << Dij << " " << Ui[0] << " " << Uj[0] << std::endl;
+ abort();
+ }
+
+ const Real UDij = Uj[0]/Dij;
+ // RR-block
+ A(i,j) = Uj[0]*UDij;
+ bnbr += Uj[1]*Uj[1]*UDij; // keep in mind: multiplication by factor 2 below!!
+ // Rp/pp/pR-block
+ const Real Uj2Dij3 = UDij*UDij/Dij;
+ // factors for matrix
+ //
+ // Rj^3
+ // ---------
+ // 2 dij^3
+ const Real RpUDij3 = static_cast<Real>(0.5)*Uj2Dij3*Uj[0];
+ //
+ // - Ri Rj^2
+ // ----------
+ // dij^3
+ const Real pRUDij3 = static_cast<Real>(-1.0)*Uj2Dij3*Ui[0];
+ //
+ // Ri Rj^3
+ // ---------
+ // 2 dij^3
+ const Real ppUDij3 = RpUDij3*Ui[0];
+ //
+ // - 3 Ri Rj^3
+ // ------------
+ // 2 dij^5
+ const Real ppUDij5 = static_cast<Real>(-3.0)*ppUDij3/(Dij*Dij);
+ // factors for rhs
+ //
+ // radius part
+ //
+ // Rj^2 Rdj
+ // --------- (comment 1/4 due to multiplication by factor 2 below)
+ // 2 dij^3
+ const Real rhsR1 = static_cast<Real>(0.25)*Uj2Dij3*Uj[1];
+ //
+ // Rj^3
+ // --------- (comment 1/8 due to multiplication by factor 2 below)
+ // 4 dij^3
+ const Real rhsR2 = static_cast<Real>(0.125)*Uj2Dij3*Uj[0];
+ //
+ // 3 Rj^3
+ // --------- (comment 3/8 due to multiplication by factor 2 below)
+ // 4 dij^5
+ const Real velpro = (Uj[2] - Ui[2])*Uj[3] + (Uj[4] - Ui[4])*Uj[5] + (Uj[6] - Ui[6])*Uj[7];
+ const Real rhsR3 = static_cast<Real>(3.0)*rhsR2*velpro/(Dij*Dij);
+ //
+ // position part
+ //
+ // Rj Rdj
+ // -------- (2 Ri Rdj + Rdi Rj)
+ // dij^3
+ const Real rhsp1 = (static_cast<Real>(2.0)*Ui[0]*Uj[1]+Ui[1]*Uj[0])*UDij*Uj[1]/(Dij*Dij);
+ //
+ // Rj^2
+ // -------- (Rdi Rj + 5 Ri Rdj)
+ // 2 dij^3
+ const Real rhsp2 = static_cast<Real>(0.5)*Uj2Dij3*(Ui[1]*Uj[0]+static_cast<Real>(5.0)*Ui[0]*Uj[1]);
+ //
+ // 3 Rj^2
+ // -------- (Rdi Rj + 5 Ri Rdj)
+ // 2 dij^5
+ const Real rhsp3 = static_cast<Real>(3.0)*rhsp2/(Dij*Dij) * velpro;
+ for (std::size_t rr=0; rr<3; rr++)
+ {
+ const int jp = _N+3*j+rr;
+ const int ip = _N+3*i+rr;
+ const Real pos_ij = Ui[rr*2+2] - Uj[rr*2+2];
+ // Rp-block
+ A(i,jp) = RpUDij3*pos_ij;
+ // pR-block
+ A(ip,j) = pRUDij3*pos_ij;
+ // pp-block
+ A(ip,jp) = ppUDij3;
+ A(ip,jp) += ppUDij5*pos_ij*pos_ij;
+ const int ia = (rr+1)%3;
+ const int ib = (rr+2)%3;
+ A(ip,_N+3*j+ia) = ppUDij5*pos_ij*(Ui[ia*2+2] - Uj[ia*2+2]);
+ A(ip,_N+3*j+ib) = ppUDij5*pos_ij*(Ui[ib*2+2] - Uj[ib*2+2]);
+
+ // rhs contributions
+ Real vel_ij = Ui[rr*2+3] + static_cast<Real>(5.0) * Uj[rr*2+3];
+ bnbr += rhsR1*pos_ij*vel_ij;
+ vel_ij = Ui[rr*2+3] + static_cast<Real>(2.0) * Uj[rr*2+3];
+ bnbr -= rhsR2*Uj[rr*2+3]*vel_ij;
+ bnbr -= rhsR3*pos_ij*vel_ij;
+ bnbrp[rr] += rhsp1*pos_ij;
+ bnbrp[rr] -= rhsp2*Uj[rr*2+3];
+ bnbrp[rr] -= rhsp3*pos_ij;
+ }
+ }
+ }
+
+ // 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];
+ // coupling terms
+ b(i) -= static_cast<Real>(2.0)*bnbr;
+ // translation term
+ b(i) += static_cast<Real>(0.25)*(Ui[3]*Ui[3] + Ui[5]*Ui[5] + Ui[7]*Ui[7]);
+
+ // pressure part
+ for (std::size_t rr=0; rr<3; rr++)
+ {
+ const int ip = _N+3*i+rr;
+ // coupling terms
+ b(ip) = bnbrp[rr];
+ // std rhs
+ b(ip) -= Ui[1]*Ui[rr*2+3];
+ // primary Bjerknes force and viscous drag neglected
+ }
+ assert(!isnan(b(i)));
+ }
+
+ //std::cout << A << std::endl << std::endl;
+ //std::cout << b << 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 State8& Ui = U[i];
+ State8& ri = rhs[i];
+ ri[0] = Ui[1];
+ ri[1] = Rddot(i);
+ ri[2] = Ui[3];
+ ri[3] = Rddot(_N+3*i);
+ ri[4] = Ui[5];
+ ri[5] = Rddot(_N+3*i+1);
+ ri[6] = Ui[7];
+ ri[7] = Rddot(_N+3*i+2);
+ }
+ }
+
+ 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;
+
+ const size_t N = 5*_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);
+ out[i+2*_N] = U[i][2];
+ out[i+3*_N] = U[i][4];
+ out[i+4*_N] = U[i][6];
+ Vg += 4.0/3.0*r*r*r*M_PI;
+ }
+ out[5*_N] = bd.pext(t);
+ out[5*_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 << U[i][2] << " " << U[i][4] << " " << U[i][6] << 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;
+
+ vtk << "VECTORS transvel double\n";
+ vtk << "LOOKUP_TABLE default" << std::endl;
+ for (int i = 0; i < bd.Nbubbles; ++i)
+ vtk << U[i][3] << " " << U[i][5] << " " << U[i][7] << std::endl;
+ vtk << std::endl;
+ }
+ }
+};
+
+#endif /* RPCLUSTERPOS_D_H_UIJM88YU */
