commit 3b14e1736fe91a2fe0b2480ee40b0ca3a8bf2c53
parent 7aa18619aa9c531bfc74550d1dbe30c34c2b22e4
Author: Ursula Rasthofer <urasthofer@ethz.ch>
Date: Tue, 7 Feb 2017 12:20:31 +0100
KM cluster with higher-order terms in 1/c according to Fuster & Colonius 2011
Diffstat:
1 file changed, 214 insertions(+), 0 deletions(-)
diff --git a/src/kernels/KMCluster_FC.h b/src/kernels/KMCluster_FC.h
@@ -0,0 +1,214 @@
+/* File: KMCluster_FC.h */
+/* Date: Wed Jan 25 15:10:28 2017 */
+/* Author: Ursula Rasthofer */
+/* Tag: Bubble cluster based on Keller-Miksis following Fuster & Colonius (2011) */
+/* Copyright 2017 ETH Zurich. All Rights Reserved. */
+#ifndef BUBBLECLUSTER_FC_H_UJH09T7V
+#define BUBBLECLUSTER_FC_H_UJH09T7V
+
+#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 KMCluster_FC : public KernelBase<Tinput,Trhs>
+{
+ const size_t _N;
+ GnuplotDump* m_dumper;
+
+public:
+ KMCluster_FC(const size_t N=0, const int format=GnuplotDump::BIN) : _N(N) { m_dumper = new GnuplotDump("out", format); }
+ virtual ~KMCluster_FC() { 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 << "KMCluster_FC: 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 << "KMCluster_FC: 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);
+
+ const Real clInv = 1.0/bd.cl;
+
+ for (size_t i = 0; i < _N; ++i)
+ {
+ const State2& Ui = U[i];
+ // fill row i of A
+ Real bnbr = 0;
+ Real bnbr2 = 0;
+ for (size_t j = 0; j < _N; ++j)
+ {
+ const State2& 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;
+ else
+ {
+ const Real UDij = Uj[0]/bd.Dij[j + i*_N];
+ A(i,j) = (static_cast<Real>(1) + Ui[1]*clInv)*Uj[0]*UDij + static_cast<Real>(6.0)*Uj[1]*UDij*Ui[0]*clInv;
+ bnbr += Uj[1]*Uj[1]*UDij;
+ bnbr2 += Uj[1]*Uj[1]*Uj[1]/bd.Dij[j + i*_N];
+ }
+ }
+
+ // 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*Rdc*bd.dpext(t + Ui[0]*clInv);
+ // b(i) -= bd.rInv*Rdc*bd.dpext(t);
+ b(i) -= static_cast<Real>(2.0)*(static_cast<Real>(1.0) + Rdc)*bnbr;
+ b(i) -= static_cast<Real>(2.0)*Ui[0]*clInv*bnbr2;
+
+ assert(!isnan(b(i)));
+ }
+
+ /* std::cout << A << std::endl << std::endl; */
+ // Eigen::Matrix<Real, Eigen::Dynamic, 1> Rddot = A.fullPivHouseholderQr().solve(b);
+ 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 /* BUBBLECLUSTER_FC_H_UJH09T7V */
