KMCluster_FC.h (7787B)
1 /* File: KMCluster_FC.h */ 2 /* Date: Wed Jan 25 15:10:28 2017 */ 3 /* Author: Ursula Rasthofer */ 4 /* Tag: Bubble cluster based on Keller-Miksis following Fuster & Colonius (2011) */ 5 /* Copyright 2017 ETH Zurich. All Rights Reserved. */ 6 #ifndef BUBBLECLUSTER_FC_H_UJH09T7V 7 #define BUBBLECLUSTER_FC_H_UJH09T7V 8 9 #include "BubbleBase.h" 10 11 #include <ODETB/GnuplotDump.h> 12 #include <ODETB/TimeStepper/KernelBase.h> 13 14 #include <Eigen/Core> 15 #include <Eigen/Dense> 16 17 #include <cassert> 18 #include <cstdlib> 19 #include <fstream> 20 #include <cmath> 21 #include <vector> 22 23 24 template <typename Tinput, typename Trhs=Tinput> 25 class KMCluster_FC : public KernelBase<Tinput,Trhs> 26 { 27 const size_t _N; 28 GnuplotDump* m_dumper; 29 30 public: 31 KMCluster_FC(const size_t N=0, const int format=GnuplotDump::BIN) : _N(N) { m_dumper = new GnuplotDump("out", format); } 32 virtual ~KMCluster_FC() { delete m_dumper;} 33 34 static void IC(Tinput& U, BubbleData& data) 35 { 36 assert(Tinput::DataType::SIZE == 2); 37 assert(data.Nbubbles == U.size()); 38 const size_t N = data.Nbubbles; 39 40 // read IC from file 41 std::ifstream infile(data.filename); 42 size_t k; 43 infile >> k; 44 if (N != k) 45 { 46 std::cerr << "KMCluster_FC: Allocated " << N << " bubbles."<< std::endl; 47 std::cerr << " File " << data.filename << " expects " << k << " bubbles."<< std::endl; 48 abort(); 49 } 50 51 for (size_t i = 0; i < N; ++i) 52 { 53 if (infile.good()) 54 { 55 Bubble& b = data.coords[i]; 56 infile >> b.pos[0] >> b.pos[1] >> b.pos[2]; 57 infile >> data.R0[i]; 58 infile >> data.Rdot0[i]; 59 60 typename Tinput::DataType& IC = U[i]; 61 IC[0] = data.R0[i]; 62 IC[1] = data.Rdot0[i]; 63 } 64 else 65 { 66 std::cerr << "KMCluster_FC: Input file has incorrect format" << std::endl; 67 abort(); 68 } 69 } 70 71 for (size_t i = 0; i < N; ++i) { 72 const Bubble& bi = data.coords[i]; 73 for (size_t j = 0; j < N; ++j) { 74 const Bubble& bj = data.coords[j]; 75 data.Dij[j + i*N] = bi.distance(bj); 76 } 77 } 78 79 #ifndef NDEBUG 80 for (size_t i = 0; i < N; ++i) 81 assert(0 == data.Dij[i + i*N]); 82 #endif /* NDEBUG */ 83 } 84 85 virtual void compute(const Tinput& U, Trhs& rhs, Real const t, void const* const data=nullptr) 86 { 87 const BubbleData& bd = *(BubbleData const* const)data; 88 89 assert(_N == U.size()); 90 Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic> A(_N,_N); 91 Eigen::Matrix<Real, Eigen::Dynamic, 1> b(_N); 92 93 const Real clInv = 1.0/bd.cl; 94 95 for (size_t i = 0; i < _N; ++i) 96 { 97 const typename Tinput::DataType& Ui = U[i]; 98 // fill row i of A 99 Real bnbr = 0; 100 Real bnbr2 = 0; 101 for (size_t j = 0; j < _N; ++j) 102 { 103 const typename Tinput::DataType& Uj = U[j]; 104 if (i == j) 105 A(i,j) = (static_cast<Real>(1) - Ui[1]*clInv)*Ui[0] + static_cast<Real>(4.0)*bd.nuL*clInv; 106 else 107 { 108 const Real UDij = Uj[0]/bd.Dij[j + i*_N]; 109 A(i,j) = (static_cast<Real>(1) + Ui[1]*clInv)*Uj[0]*UDij + static_cast<Real>(1.0)*Uj[1]*UDij*Ui[0]*clInv; 110 // extension of Matrix: keeping all terms would yield a factor of 6 for the second term, 111 // following simplifications of Fuster and Colonius results in a factor of 1 112 bnbr += Uj[1]*Uj[1]*UDij; 113 bnbr2 += Uj[1]*Uj[1]*Uj[1]/bd.Dij[j + i*_N]; 114 } 115 } 116 117 // compute bi 118 const Real Rdc = Ui[1]*clInv; 119 const Real pB = (bd.p0 - bd.pv + static_cast<Real>(2.0)*bd.sigma/bd.R0[i])*std::pow(bd.R0[i]/Ui[0], 3.0*bd.gamma); // bubble pressure 120 121 b(i) = (static_cast<Real>(0.5)*Rdc - static_cast<Real>(1.5))*Ui[1]*Ui[1]; 122 b(i) += bd.rInv*(static_cast<Real>(1.0) + static_cast<Real>(1.0 - 3.0*bd.gamma)*Rdc)*pB; 123 b(i) -= bd.rInv*static_cast<Real>(2.0)*bd.sigma/Ui[0]; 124 b(i) -= static_cast<Real>(4.0)*bd.nuL*Ui[1]/Ui[0]; 125 b(i) -= bd.rInv*(static_cast<Real>(1.0) + Rdc)*(bd.p0 - bd.pv + bd.pext(t + Ui[0]*clInv)); 126 // b(i) -= bd.rInv*(static_cast<Real>(1.0) + Rdc)*(bd.p0 - bd.pv + bd.pext(t)); 127 b(i) -= bd.rInv*Ui[0]*clInv*bd.dpext(t + Ui[0]*clInv); 128 // b(i) -= bd.rInv*Ui[0]*clInv*bd.dpext(t); 129 b(i) -= static_cast<Real>(2.0)*(static_cast<Real>(1.0) + Rdc)*bnbr; 130 b(i) -= static_cast<Real>(2.0)*Ui[0]*clInv*bnbr2; 131 132 assert(!isnan(b(i))); 133 } 134 135 /* std::cout << A << std::endl << std::endl; */ 136 // Eigen::Matrix<Real, Eigen::Dynamic, 1> Rddot = A.fullPivHouseholderQr().solve(b); 137 Eigen::Matrix<Real, Eigen::Dynamic, 1> Rddot = A.colPivHouseholderQr().solve(b); 138 139 for (size_t i = 0; i < _N; ++i) 140 { 141 const typename Tinput::DataType& Ui = U[i]; 142 typename Trhs::DataType& ri = rhs[i]; 143 ri[0] = Ui[1]; 144 ri[1] = Rddot(i); 145 } 146 } 147 148 virtual void write(const size_t step, const Real t, const Real dt, const Tinput& U, const void * const data=nullptr) 149 { 150 static size_t callCount = 0; 151 152 const BubbleData& bd = *(BubbleData const* const)data; 153 154 // double pb[3]; 155 // const int id[3] = {0, 7, 1}; 156 // for (int i = 0; i < 3; ++i) 157 // pb[i] = bd.pBubble(U, id[i]); 158 159 // // class 1, 2, 3, 4, 5, 6 160 // // const Real odat[6] = {U[0][0], U[1][0], U[4][0], U[2][0], U[5][0], U[3][0]}; 161 // const Real odat[7] = {U[0][0], U[7][0], U[1][0], pb[0], pb[1], pb[2], bd.pext(t)}; 162 163 const size_t N = 2*_N+2; 164 std::vector<double> out(N); 165 double Vg = 0.0; 166 for (size_t i = 0; i < _N; ++i) { 167 const double r = U[i][0]; 168 out[i] = r; 169 out[i+_N] = bd.pBubble(U, i); 170 Vg += 4.0/3.0*r*r*r*M_PI; 171 } 172 out[2*_N] = bd.pext(t); 173 out[2*_N+1] = Vg; 174 175 m_dumper->write(step, t, dt, out.data(), N); 176 177 if (bd.bVTK) 178 { 179 // dump bubble vtk 180 std::stringstream streamer; 181 streamer<< "bubbles-vtk_" << setfill('0') << setw(8) << callCount++ << ".vtk"; 182 std::ofstream vtk(streamer.str().c_str()); 183 vtk.setf(std::ios::scientific, std::ios::floatfield); 184 185 vtk << "# vtk DataFile Version 2.0\n"; 186 vtk << "Unstructured grid legacy vtk file with point scalar data\n"; 187 vtk << "ASCII\n" << std::endl; 188 189 vtk << "DATASET UNSTRUCTURED_GRID\n"; 190 vtk << "POINTS " << bd.Nbubbles << " double" << std::endl; 191 for (size_t i = 0; i < bd.Nbubbles; ++i) 192 vtk << bd.coords[i].pos[0] << " " << bd.coords[i].pos[1] << " " << bd.coords[i].pos[2] << std::endl; 193 vtk << std::endl; 194 195 vtk << "POINT_DATA " << bd.Nbubbles << std::endl; 196 vtk << "SCALARS radius double\n"; 197 vtk << "LOOKUP_TABLE default" << std::endl; 198 for (size_t i = 0; i < bd.Nbubbles; ++i) 199 vtk << U[i][0] << std::endl; 200 vtk << std::endl; 201 202 // vtk << "POINT_DATA " << bd.Nbubbles << std::endl; 203 vtk << "SCALARS velocity double\n"; 204 vtk << "LOOKUP_TABLE default" << std::endl; 205 for (size_t i = 0; i < bd.Nbubbles; ++i) 206 vtk << U[i][1] << std::endl; 207 vtk << std::endl; 208 } 209 } 210 }; 211 212 #endif /* BUBBLECLUSTER_FC_H_UJH09T7V */