RPClusterPositions_D.h (13016B)
1 /* File: RPClusterPositions_D.h */ 2 /* Date: Tue Feb 7 09:26:31 2017 */ 3 /* Author: Ursula Rasthofer */ 4 /* Tag: Bubble cluster based on Rayleigh-Plesset */ 5 /* including translation of bubbles following Doinikov (2004) */ 6 /* Copyright 2017 ETH Zurich. All Rights Reserved. */ 7 #ifndef RPCLUSTERPOS_D_H_UIJM88YU 8 #define RPCLUSTERPOS_D_H_UIJM88YU 9 10 #include "BubbleBase.h" 11 12 #include <ODETB/GnuplotDump.h> 13 #include <ODETB/TimeStepper/KernelBase.h> 14 15 #include <Eigen/Core> 16 #include <Eigen/Dense> 17 18 #include <cassert> 19 #include <cstdlib> 20 #include <fstream> 21 #include <cmath> 22 #include <vector> 23 24 25 template <typename Tinput, typename Trhs=Tinput> 26 class RPClusterPositions_D : public KernelBase<Tinput,Trhs> 27 { 28 const size_t _N; 29 GnuplotDump* m_dumper; 30 31 public: 32 RPClusterPositions_D(const size_t N=0, const int format=GnuplotDump::BIN) : _N(N) { m_dumper = new GnuplotDump("out", format); } 33 virtual ~RPClusterPositions_D() { delete m_dumper;} 34 35 static void IC(Tinput& U, BubbleData& data) 36 { 37 assert(Tinput::DataType::SIZE == 8); 38 assert(data.Nbubbles == U.size()); 39 const size_t N = data.Nbubbles; 40 41 // read IC from file 42 std::ifstream infile(data.filename); 43 size_t k; 44 infile >> k; 45 if (N != k) 46 { 47 std::cerr << "RPClusterPositions_D: Allocated " << N << " bubbles."<< std::endl; 48 std::cerr << " File " << data.filename << " expects " << k << " bubbles."<< std::endl; 49 abort(); 50 } 51 52 for (size_t i = 0; i < N; ++i) 53 { 54 if (infile.good()) 55 { 56 Bubble& b = data.coords[i]; 57 infile >> b.pos[0] >> b.pos[1] >> b.pos[2]; 58 infile >> data.R0[i]; 59 infile >> data.Rdot0[i]; 60 61 typename Tinput::DataType& IC = U[i]; 62 IC[0] = data.R0[i]; 63 IC[1] = data.Rdot0[i]; 64 IC[2] = b.pos[0]; 65 IC[3] = 0.0; 66 IC[4] = b.pos[1]; 67 IC[5] = 0.0; 68 IC[6] = b.pos[2]; 69 IC[7] = 0.0; 70 } 71 else 72 { 73 std::cerr << "RPClusterPositions_D: Input file has incorrect format" << std::endl; 74 abort(); 75 } 76 } 77 78 for (size_t i = 0; i < N; ++i) { 79 const Bubble& bi = data.coords[i]; 80 for (size_t j = 0; j < N; ++j) { 81 const Bubble& bj = data.coords[j]; 82 data.Dij[j + i*N] = bi.distance(bj); 83 } 84 } 85 86 #ifndef NDEBUG 87 for (size_t i = 0; i < N; ++i) 88 assert(0 == data.Dij[i + i*N]); 89 #endif /* NDEBUG */ 90 } 91 92 virtual void compute(const Tinput& U, Trhs& rhs, Real const t, void const* const data=nullptr) 93 { 94 const BubbleData& bd = *(BubbleData const* const)data; 95 96 assert(_N == U.size()); 97 const int size = 4*_N; 98 Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic> A(size,size); 99 Eigen::Matrix<Real, Eigen::Dynamic, 1> b(size); 100 101 for (size_t i = 0; i < _N; ++i) 102 { 103 const typename Tinput::DataType& Ui = U[i]; 104 const Bubble bi = {Ui[2], Ui[4], Ui[6]}; 105 // fill row i of A 106 Real bnbr = 0; 107 std::vector<Real> bnbrp(3,0.0); 108 for (size_t j = 0; j < _N; ++j) 109 { 110 const typename Tinput::DataType& Uj = U[j]; 111 if (i == j) 112 { 113 A(i,j) = Ui[0]; 114 for (std::size_t rr=0; rr<3; rr++) 115 { 116 const int ip = _N+3*i+rr; 117 A(i,ip) = static_cast<Real>(0.0); 118 A(ip,ip) = Ui[0]/static_cast<Real>(3.0); 119 A(ip,_N+3*i+(rr+1)%3) = static_cast<Real>(0.0); 120 A(ip,_N+3*i+(rr+2)%3) = static_cast<Real>(0.0); 121 A(ip,j) = static_cast<Real>(0.0); 122 } 123 } 124 else 125 { 126 // compute distance between bubbles 127 const Bubble bj = {Uj[2], Uj[4], Uj[6]}; 128 const Real Dij = bi.distance(bj); 129 130 if (Dij < (Ui[0]+Uj[0])) 131 { 132 std::cerr << "Colliding bubbles " << Dij << " " << Ui[0] << " " << Uj[0] << std::endl; 133 abort(); 134 } 135 136 const Real UDij = Uj[0]/Dij; 137 // RR-block 138 A(i,j) = Uj[0]*UDij; 139 bnbr += Uj[1]*Uj[1]*UDij; // keep in mind: multiplication by factor 2 below!! 140 // Rp/pp/pR-block 141 const Real Uj2Dij3 = UDij*UDij/Dij; 142 // factors for matrix 143 // 144 // Rj^3 145 // --------- 146 // 2 dij^3 147 const Real RpUDij3 = static_cast<Real>(0.5)*Uj2Dij3*Uj[0]; 148 // 149 // - Ri Rj^2 150 // ---------- 151 // dij^3 152 const Real pRUDij3 = static_cast<Real>(-1.0)*Uj2Dij3*Ui[0]; 153 // 154 // Ri Rj^3 155 // --------- 156 // 2 dij^3 157 const Real ppUDij3 = RpUDij3*Ui[0]; 158 // 159 // - 3 Ri Rj^3 160 // ------------ 161 // 2 dij^5 162 const Real ppUDij5 = static_cast<Real>(-3.0)*ppUDij3/(Dij*Dij); 163 // factors for rhs 164 // 165 // radius part 166 // 167 // Rj^2 Rdj 168 // --------- (comment 1/4 due to multiplication by factor 2 below) 169 // 2 dij^3 170 const Real rhsR1 = static_cast<Real>(0.25)*Uj2Dij3*Uj[1]; 171 // 172 // Rj^3 173 // --------- (comment 1/8 due to multiplication by factor 2 below) 174 // 4 dij^3 175 const Real rhsR2 = static_cast<Real>(0.125)*Uj2Dij3*Uj[0]; 176 // 177 // 3 Rj^3 178 // --------- (comment 3/8 due to multiplication by factor 2 below) 179 // 4 dij^5 180 const Real velpro = (Uj[2] - Ui[2])*Uj[3] + (Uj[4] - Ui[4])*Uj[5] + (Uj[6] - Ui[6])*Uj[7]; 181 const Real rhsR3 = static_cast<Real>(3.0)*rhsR2*velpro/(Dij*Dij); 182 // 183 // position part 184 // 185 // Rj Rdj 186 // -------- (2 Ri Rdj + Rdi Rj) 187 // dij^3 188 const Real rhsp1 = (static_cast<Real>(2.0)*Ui[0]*Uj[1]+Ui[1]*Uj[0])*UDij*Uj[1]/(Dij*Dij); 189 // 190 // Rj^2 191 // -------- (Rdi Rj + 5 Ri Rdj) 192 // 2 dij^3 193 const Real rhsp2 = static_cast<Real>(0.5)*Uj2Dij3*(Ui[1]*Uj[0]+static_cast<Real>(5.0)*Ui[0]*Uj[1]); 194 // 195 // 3 Rj^2 196 // -------- (Rdi Rj + 5 Ri Rdj) 197 // 2 dij^5 198 const Real rhsp3 = static_cast<Real>(3.0)*rhsp2/(Dij*Dij) * velpro; 199 for (std::size_t rr=0; rr<3; rr++) 200 { 201 const int jp = _N+3*j+rr; 202 const int ip = _N+3*i+rr; 203 const Real pos_ij = Ui[rr*2+2] - Uj[rr*2+2]; 204 // Rp-block 205 A(i,jp) = RpUDij3*pos_ij; 206 // pR-block 207 A(ip,j) = pRUDij3*pos_ij; 208 // pp-block 209 A(ip,jp) = ppUDij3; 210 A(ip,jp) += ppUDij5*pos_ij*pos_ij; 211 const int ia = (rr+1)%3; 212 const int ib = (rr+2)%3; 213 A(ip,_N+3*j+ia) = ppUDij5*pos_ij*(Ui[ia*2+2] - Uj[ia*2+2]); 214 A(ip,_N+3*j+ib) = ppUDij5*pos_ij*(Ui[ib*2+2] - Uj[ib*2+2]); 215 216 // rhs contributions 217 Real vel_ij = Ui[rr*2+3] + static_cast<Real>(5.0) * Uj[rr*2+3]; 218 bnbr += rhsR1*pos_ij*vel_ij; 219 vel_ij = Ui[rr*2+3] + static_cast<Real>(2.0) * Uj[rr*2+3]; 220 bnbr -= rhsR2*Uj[rr*2+3]*vel_ij; 221 bnbr -= rhsR3*pos_ij*vel_ij; 222 bnbrp[rr] += rhsp1*pos_ij; 223 bnbrp[rr] -= rhsp2*Uj[rr*2+3]; 224 bnbrp[rr] -= rhsp3*pos_ij; 225 } 226 } 227 } 228 229 // compute bi 230 const Real pstat = bd.p0 - bd.pv; 231 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); 232 233 b(i) = 0.0; 234 b(i) -= static_cast<Real>(1.5)*Ui[1]*Ui[1]; 235 b(i) -= bd.rInv*(pstat + bd.pext(t)); 236 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]; 237 // coupling terms 238 b(i) -= static_cast<Real>(2.0)*bnbr; 239 // translation term 240 b(i) += static_cast<Real>(0.25)*(Ui[3]*Ui[3] + Ui[5]*Ui[5] + Ui[7]*Ui[7]); 241 242 // pressure part 243 for (std::size_t rr=0; rr<3; rr++) 244 { 245 const int ip = _N+3*i+rr; 246 // coupling terms 247 b(ip) = bnbrp[rr]; 248 // std rhs 249 b(ip) -= Ui[1]*Ui[rr*2+3]; 250 // primary Bjerknes force and viscous drag neglected 251 } 252 assert(!isnan(b(i))); 253 } 254 255 //std::cout << A << std::endl << std::endl; 256 //std::cout << b << std::endl << std::endl; 257 //Eigen::Matrix<Real, Eigen::Dynamic, 1> Rddot = A.fullPivHouseholderQr().solve(b); 258 // this is notably faster and sufficiently accurate (rasthofer Feb 3, 2017) 259 Eigen::Matrix<Real, Eigen::Dynamic, 1> Rddot = A.colPivHouseholderQr().solve(b); 260 261 for (size_t i = 0; i < _N; ++i) 262 { 263 const typename Tinput::DataType& Ui = U[i]; 264 typename Trhs::DataType& ri = rhs[i]; 265 ri[0] = Ui[1]; 266 ri[1] = Rddot(i); 267 ri[2] = Ui[3]; 268 ri[3] = Rddot(_N+3*i); 269 ri[4] = Ui[5]; 270 ri[5] = Rddot(_N+3*i+1); 271 ri[6] = Ui[7]; 272 ri[7] = Rddot(_N+3*i+2); 273 } 274 } 275 276 virtual void write(const size_t step, const Real t, const Real dt, const Tinput& U, const void * const data=nullptr) 277 { 278 static size_t callCount = 0; 279 280 const BubbleData& bd = *(BubbleData const* const)data; 281 282 const size_t N = 5*_N+2; 283 std::vector<double> out(N); 284 double Vg = 0.0; 285 for (size_t i = 0; i < _N; ++i) { 286 const double r = U[i][0]; 287 out[i] = r; 288 out[i+_N] = bd.pBubble(U, i); 289 out[i+2*_N] = U[i][2]; 290 out[i+3*_N] = U[i][4]; 291 out[i+4*_N] = U[i][6]; 292 Vg += 4.0/3.0*r*r*r*M_PI; 293 } 294 out[5*_N] = bd.pext(t); 295 out[5*_N+1] = Vg; 296 297 m_dumper->write(step, t, dt, out.data(), N); 298 299 if (bd.bVTK) 300 { 301 // dump bubble vtk 302 std::stringstream streamer; 303 streamer<< "bubbles-vtk_" << setfill('0') << setw(8) << callCount++ << ".vtk"; 304 std::ofstream vtk(streamer.str().c_str()); 305 vtk.setf(std::ios::scientific, std::ios::floatfield); 306 307 vtk << "# vtk DataFile Version 2.0\n"; 308 vtk << "Unstructured grid legacy vtk file with point scalar data\n"; 309 vtk << "ASCII\n" << std::endl; 310 311 vtk << "DATASET UNSTRUCTURED_GRID\n"; 312 vtk << "POINTS " << bd.Nbubbles << " double" << std::endl; 313 for (size_t i = 0; i < bd.Nbubbles; ++i) 314 vtk << U[i][2] << " " << U[i][4] << " " << U[i][6] << std::endl; 315 vtk << std::endl; 316 317 vtk << "POINT_DATA " << bd.Nbubbles << std::endl; 318 vtk << "SCALARS radius double\n"; 319 vtk << "LOOKUP_TABLE default" << std::endl; 320 for (size_t i = 0; i < bd.Nbubbles; ++i) 321 vtk << U[i][0] << std::endl; 322 vtk << std::endl; 323 324 // vtk << "POINT_DATA " << bd.Nbubbles << std::endl; 325 vtk << "SCALARS velocity double\n"; 326 vtk << "LOOKUP_TABLE default" << std::endl; 327 for (size_t i = 0; i < bd.Nbubbles; ++i) 328 vtk << U[i][1] << std::endl; 329 vtk << std::endl; 330 331 vtk << "VECTORS transvel double\n"; 332 vtk << "LOOKUP_TABLE default" << std::endl; 333 for (size_t i = 0; i < bd.Nbubbles; ++i) 334 vtk << U[i][3] << " " << U[i][5] << " " << U[i][7] << std::endl; 335 vtk << std::endl; 336 } 337 } 338 }; 339 340 #endif /* RPCLUSTERPOS_D_H_UIJM88YU */