RPCluster.h (6799B)
1 /* File: RPCluster.h */ 2 /* Date: Sun Feb 05 22:10:34 2017 */ 3 /* Author: Ursula Rasthofer */ 4 /* Tag: Bubble cluster based on Rayleigh-Plesset */ 5 /* Copyright 2017 ETH Zurich. All Rights Reserved. */ 6 #ifndef RPCLUSTER_H_67YUOP08 7 #define RPCLUSTER_H_67YUOP08 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 RPCluster : public KernelBase<Tinput,Trhs> 26 { 27 const size_t _N; 28 GnuplotDump* m_dumper; 29 30 public: 31 RPCluster(const size_t N=0, const int format=GnuplotDump::BIN) : _N(N) { m_dumper = new GnuplotDump("out", format); } 32 virtual ~RPCluster() { 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 << "RPCluster: 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 << "RPCluster: 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 for (size_t i = 0; i < _N; ++i) 94 { 95 const typename Tinput::DataType& Ui = U[i]; 96 // fill row i of A 97 Real bnbr = 0; 98 for (size_t j = 0; j < _N; ++j) 99 { 100 const typename Tinput::DataType& Uj = U[j]; 101 if (i == j) 102 A(i,j) = Ui[0]; 103 else 104 { 105 const Real UDij = Uj[0]/bd.Dij[j + i*_N]; 106 A(i,j) = Uj[0]*UDij; 107 bnbr += Uj[1]*Uj[1]*UDij; 108 } 109 } 110 111 // compute bi 112 const Real pstat = bd.p0 - bd.pv; 113 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); 114 115 b(i) = 0.0; 116 b(i) -= static_cast<Real>(1.5)*Ui[1]*Ui[1]; 117 b(i) -= bd.rInv*(pstat + bd.pext(t)); 118 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]; 119 b(i) -= static_cast<Real>(2.0)*bnbr; 120 121 assert(!isnan(b(i))); 122 } 123 124 /* std::cout << A << std::endl << std::endl; */ 125 //Eigen::Matrix<Real, Eigen::Dynamic, 1> Rddot = A.fullPivHouseholderQr().solve(b); 126 // this is notably faster and sufficiently accurate (rasthofer Feb 3, 2017) 127 Eigen::Matrix<Real, Eigen::Dynamic, 1> Rddot = A.colPivHouseholderQr().solve(b); 128 129 for (size_t i = 0; i < _N; ++i) 130 { 131 const typename Tinput::DataType& Ui = U[i]; 132 typename Trhs::DataType& ri = rhs[i]; 133 ri[0] = Ui[1]; 134 ri[1] = Rddot(i); 135 } 136 } 137 138 virtual void write(const size_t step, const Real t, const Real dt, const Tinput& U, const void * const data=nullptr) 139 { 140 static size_t callCount = 0; 141 142 const BubbleData& bd = *(BubbleData const* const)data; 143 144 // double pb[3]; 145 // const int id[3] = {0, 7, 1}; 146 // for (int i = 0; i < 3; ++i) 147 // pb[i] = bd.pBubble(U, id[i]); 148 149 // // class 1, 2, 3, 4, 5, 6 150 // // const Real odat[6] = {U[0][0], U[1][0], U[4][0], U[2][0], U[5][0], U[3][0]}; 151 // const Real odat[7] = {U[0][0], U[7][0], U[1][0], pb[0], pb[1], pb[2], bd.pext(t)}; 152 153 const size_t N = 2*_N+2; 154 std::vector<double> out(N); 155 double Vg = 0.0; 156 for (size_t i = 0; i < _N; ++i) { 157 const double r = U[i][0]; 158 out[i] = r; 159 out[i+_N] = bd.pBubble(U, i); 160 Vg += 4.0/3.0*r*r*r*M_PI; 161 } 162 out[2*_N] = bd.pext(t); 163 out[2*_N+1] = Vg; 164 165 m_dumper->write(step, t, dt, out.data(), N); 166 167 if (bd.bVTK) 168 { 169 // dump bubble vtk 170 std::stringstream streamer; 171 streamer<< "bubbles-vtk_" << setfill('0') << setw(8) << callCount++ << ".vtk"; 172 std::ofstream vtk(streamer.str().c_str()); 173 vtk.setf(std::ios::scientific, std::ios::floatfield); 174 175 vtk << "# vtk DataFile Version 2.0\n"; 176 vtk << "Unstructured grid legacy vtk file with point scalar data\n"; 177 vtk << "ASCII\n" << std::endl; 178 179 vtk << "DATASET UNSTRUCTURED_GRID\n"; 180 vtk << "POINTS " << bd.Nbubbles << " double" << std::endl; 181 for (size_t i = 0; i < bd.Nbubbles; ++i) 182 vtk << bd.coords[i].pos[0] << " " << bd.coords[i].pos[1] << " " << bd.coords[i].pos[2] << std::endl; 183 vtk << std::endl; 184 185 vtk << "POINT_DATA " << bd.Nbubbles << std::endl; 186 vtk << "SCALARS radius double\n"; 187 vtk << "LOOKUP_TABLE default" << std::endl; 188 for (size_t i = 0; i < bd.Nbubbles; ++i) 189 vtk << U[i][0] << std::endl; 190 vtk << std::endl; 191 192 // vtk << "POINT_DATA " << bd.Nbubbles << std::endl; 193 vtk << "SCALARS velocity double\n"; 194 vtk << "LOOKUP_TABLE default" << std::endl; 195 for (size_t i = 0; i < bd.Nbubbles; ++i) 196 vtk << U[i][1] << std::endl; 197 vtk << std::endl; 198 } 199 } 200 }; 201 202 #endif /* RPCLUSTER_H_67YUOP08 */