KMClusterPositions_D.h (14307B)
1 /* File: KMClusterPositions_D.h */ 2 /* Date: Sun Feb 5 13:24:54 2017 */ 3 /* Author: Ursula Rasthofer */ 4 /* Tag: Bubble cluster based on Keller-Miksis */ 5 /* including translation of bubbles following Doinikov (2004) */ 6 /* Copyright 2017 ETH Zurich. All Rights Reserved. */ 7 #ifndef BUBBLECLUSTERPOS_D_H_GHBNU569 8 #define BUBBLECLUSTERPOS_D_H_GHBNU569 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 <iostream> 21 #include <fstream> 22 #include <cmath> 23 #include <vector> 24 25 template <typename Tinput, typename Trhs=Tinput> 26 class KMClusterPositions_D : public KernelBase<Tinput,Trhs> 27 { 28 const size_t _N; 29 GnuplotDump* m_dumper; 30 31 public: 32 KMClusterPositions_D(const size_t N=0, const int format=GnuplotDump::BIN) : _N(N) { m_dumper = new GnuplotDump("out", format); } 33 virtual ~KMClusterPositions_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 << "KMClusterPositions_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 << "KMClusterPositions_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 const Real clInv = 1.0/bd.cl; 102 103 for (size_t i = 0; i < _N; ++i) 104 { 105 const typename Tinput::DataType& Ui = U[i]; 106 const Bubble bi = {Ui[2], Ui[4], Ui[6]}; 107 // fill row i of A 108 Real bnbr = 0; 109 std::vector<Real> bnbrp(3,0.0); 110 for (size_t j = 0; j < _N; ++j) 111 { 112 const typename Tinput::DataType& Uj = U[j]; 113 if (i == j) 114 { 115 A(i,j) = (static_cast<Real>(1) - Ui[1]*clInv)*Ui[0] + static_cast<Real>(4.0)*bd.nuL*clInv; 116 for (std::size_t rr=0; rr<3; rr++) 117 { 118 const int ip = _N+3*i+rr; 119 A(i,ip) = static_cast<Real>(0.0); 120 A(ip,ip) = Ui[0]/static_cast<Real>(3.0); 121 A(ip,_N+3*i+(rr+1)%3) = static_cast<Real>(0.0); 122 A(ip,_N+3*i+(rr+2)%3) = static_cast<Real>(0.0); 123 A(ip,j) = static_cast<Real>(0.0); 124 } 125 } 126 else 127 { 128 // compute distance between bubbles 129 const Bubble bj = {Uj[2], Uj[4], Uj[6]}; 130 const Real Dij = bi.distance(bj); 131 132 if (Dij < (Ui[0]+Uj[0])) 133 { 134 std::cerr << "Colliding bubbles " << i << " " << j << " distance " << Dij << " radius i " << Ui[0] << " radius j " << Uj[0] << std::endl; 135 abort(); 136 } 137 138 const Real UDij = Uj[0]/Dij; 139 // RR-block 140 A(i,j) = Uj[0]*UDij; 141 bnbr += Uj[1]*Uj[1]*UDij; // keep in mind: multiplication by factor 2 below!! 142 // Rp/pp/pR-block 143 const Real Uj2Dij3 = UDij*UDij/Dij; 144 // factors for matrix 145 // 146 // Rj^3 147 // --------- 148 // 2 dij^3 149 const Real RpUDij3 = static_cast<Real>(0.5)*Uj2Dij3*Uj[0]; 150 // 151 // - Ri Rj^2 152 // ---------- 153 // dij^3 154 const Real pRUDij3 = static_cast<Real>(-1.0)*Uj2Dij3*Ui[0]; 155 // 156 // Ri Rj^3 157 // --------- 158 // 2 dij^3 159 const Real ppUDij3 = RpUDij3*Ui[0]; 160 // 161 // - 3 Ri Rj^3 162 // ------------ 163 // 2 dij^5 164 const Real ppUDij5 = static_cast<Real>(-3.0)*ppUDij3/(Dij*Dij); 165 // factors for rhs 166 // 167 // radius part 168 // 169 // Rj^2 Rdj 170 // --------- (comment 1/4 due to multiplication by factor 2 below) 171 // 2 dij^3 172 const Real rhsR1 = static_cast<Real>(0.25)*Uj2Dij3*Uj[1]; 173 // 174 // Rj^3 175 // --------- (comment 1/8 due to multiplication by factor 2 below) 176 // 4 dij^3 177 const Real rhsR2 = static_cast<Real>(0.125)*Uj2Dij3*Uj[0]; 178 // 179 // 3 Rj^3 180 // --------- (comment 3/8 due to multiplication by factor 2 below) 181 // 4 dij^5 182 const Real velpro = (Uj[2] - Ui[2])*Uj[3] + (Uj[4] - Ui[4])*Uj[5] + (Uj[6] - Ui[6])*Uj[7]; 183 const Real rhsR3 = static_cast<Real>(3.0)*rhsR2*velpro/(Dij*Dij); 184 // 185 // position part 186 // 187 // Rj Rdj 188 // -------- (2 Ri Rdj + Rdi Rj) 189 // dij^3 190 const Real rhsp1 = (static_cast<Real>(2.0)*Ui[0]*Uj[1]+Ui[1]*Uj[0])*UDij*Uj[1]/(Dij*Dij); 191 // 192 // Rj^2 193 // -------- (Rdi Rj + 5 Ri Rdj) 194 // 2 dij^3 195 const Real rhsp2 = static_cast<Real>(0.5)*Uj2Dij3*(Ui[1]*Uj[0]+static_cast<Real>(5.0)*Ui[0]*Uj[1]); 196 // 197 // 3 Rj^2 198 // -------- (Rdi Rj + 5 Ri Rdj) 199 // 2 dij^5 200 const Real rhsp3 = static_cast<Real>(3.0)*rhsp2/(Dij*Dij) * velpro; // consider opposite sign of velpro below (substraction instead of addition) 201 for (std::size_t rr=0; rr<3; rr++) 202 { 203 const int jp = _N+3*j+rr; 204 const int ip = _N+3*i+rr; 205 const Real pos_ij = Ui[rr*2+2] - Uj[rr*2+2]; 206 // Rp-block 207 A(i,jp) = RpUDij3*pos_ij; 208 // pR-block 209 A(ip,j) = pRUDij3*pos_ij; 210 // pp-block 211 A(ip,jp) = ppUDij3; 212 A(ip,jp) += ppUDij5*pos_ij*pos_ij; 213 const int ia = (rr+1)%3; 214 const int ib = (rr+2)%3; 215 A(ip,_N+3*j+ia) = ppUDij5*pos_ij*(Ui[ia*2+2] - Uj[ia*2+2]); 216 A(ip,_N+3*j+ib) = ppUDij5*pos_ij*(Ui[ib*2+2] - Uj[ib*2+2]); 217 218 // rhs contributions 219 Real vel_ij = Ui[rr*2+3] + static_cast<Real>(5.0) * Uj[rr*2+3]; 220 bnbr += rhsR1*pos_ij*vel_ij; 221 vel_ij = Ui[rr*2+3] + static_cast<Real>(2.0) * Uj[rr*2+3]; 222 bnbr -= rhsR2*Uj[rr*2+3]*vel_ij; 223 bnbr -= rhsR3*pos_ij*vel_ij; 224 bnbrp[rr] += rhsp1*pos_ij; 225 bnbrp[rr] -= rhsp2*Uj[rr*2+3]; 226 bnbrp[rr] -= rhsp3*pos_ij; //TODO: save multiplication here precomputing (rhsp1-rhsp3) 227 } 228 } 229 } 230 231 // compute bi 232 const Real Rdc = Ui[1]*clInv; 233 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 234 235 b(i) = (static_cast<Real>(0.5)*Rdc - static_cast<Real>(1.5))*Ui[1]*Ui[1]; 236 b(i) += bd.rInv*(static_cast<Real>(1.0) + static_cast<Real>(1.0 - 3.0*bd.gamma)*Rdc)*pB; 237 b(i) -= bd.rInv*static_cast<Real>(2.0)*bd.sigma/Ui[0]; 238 b(i) -= static_cast<Real>(4.0)*bd.nuL*Ui[1]/Ui[0]; 239 b(i) -= bd.rInv*(static_cast<Real>(1.0) + Rdc)*(bd.p0 - bd.pv + bd.pext(t + Ui[0]*clInv)); 240 // b(i) -= bd.rInv*(static_cast<Real>(1.0) + Rdc)*(bd.p0 - bd.pv + bd.pext(t)); 241 b(i) -= bd.rInv*Ui[0]*clInv*bd.dpext(t + Ui[0]*clInv); 242 // b(i) -= bd.rInv*Ui[0]*clInv*bd.dpext(t); 243 244 // coupling terms 245 b(i) -= static_cast<Real>(2.0)*bnbr; 246 // translation term 247 b(i) += static_cast<Real>(0.25)*(Ui[3]*Ui[3] + Ui[5]*Ui[5] + Ui[7]*Ui[7]); 248 249 // pressure part 250 for (std::size_t rr=0; rr<3; rr++) 251 { 252 const int ip = _N+3*i+rr; 253 // coupling terms 254 b(ip) = bnbrp[rr]; 255 // std rhs 256 b(ip) -= Ui[1]*Ui[rr*2+3]; 257 // primary Bjerknes force and viscous drag neglected 258 } 259 assert(!isnan(b(i))); 260 } 261 262 //std::cout << A << std::endl << std::endl; 263 //std::cout << b << std::endl << std::endl; 264 //Eigen::Matrix<Real, Eigen::Dynamic, 1> Rddot = A.fullPivHouseholderQr().solve(b); 265 // this is notably faster and sufficiently accurate (rasthofer Feb 3, 2017) 266 Eigen::Matrix<Real, Eigen::Dynamic, 1> Rddot = A.colPivHouseholderQr().solve(b); 267 // don't use these solvers although they are fast (err e-2) 268 // to increase speed use less expensive time integration scheme 269 //Eigen::Matrix<Real, Eigen::Dynamic, 1> Rddot = A.llt().solve(b); 270 //Eigen::Matrix<Real, Eigen::Dynamic, 1> Rddot = A.ldlt().solve(b); 271 #ifndef NDEBUG 272 Real err = (A*Rddot - b).norm() / b.norm(); 273 274 ofstream errfile; 275 errfile.open("error.dat", ios::out | ios::app); 276 errfile << scientific << setprecision(6) << err << endl; 277 errfile.close(); 278 #endif /* NDEBUG */ 279 280 for (size_t i = 0; i < _N; ++i) 281 { 282 const typename Tinput::DataType& Ui = U[i]; 283 typename Trhs::DataType& ri = rhs[i]; 284 ri[0] = Ui[1]; 285 ri[1] = Rddot(i); 286 ri[2] = Ui[3]; 287 ri[3] = Rddot(_N+3*i); 288 ri[4] = Ui[5]; 289 ri[5] = Rddot(_N+3*i+1); 290 ri[6] = Ui[7]; 291 ri[7] = Rddot(_N+3*i+2); 292 } 293 } 294 295 virtual void write(const size_t step, const Real t, const Real dt, const Tinput& U, const void * const data=nullptr) 296 { 297 static size_t callCount = 0; 298 299 const BubbleData& bd = *(BubbleData const* const)data; 300 301 const size_t N = 5*_N+2; 302 std::vector<double> out(N); 303 double Vg = 0.0; 304 for (size_t i = 0; i < _N; ++i) { 305 const double r = U[i][0]; 306 out[i] = r; 307 out[i+_N] = bd.pBubble(U, i); 308 out[i+2*_N] = U[i][2]; 309 out[i+3*_N] = U[i][4]; 310 out[i+4*_N] = U[i][6]; 311 Vg += 4.0/3.0*r*r*r*M_PI; 312 } 313 out[5*_N] = bd.pext(t); 314 out[5*_N+1] = Vg; 315 316 m_dumper->write(step, t, dt, out.data(), N); 317 318 if (bd.bVTK) 319 { 320 // dump bubble vtk 321 std::stringstream streamer; 322 streamer<< "bubbles-vtk_" << setfill('0') << setw(8) << callCount++ << ".vtk"; 323 std::ofstream vtk(streamer.str().c_str()); 324 vtk.setf(std::ios::scientific, std::ios::floatfield); 325 326 vtk << "# vtk DataFile Version 2.0\n"; 327 vtk << "Unstructured grid legacy vtk file with point scalar data\n"; 328 vtk << "ASCII\n" << std::endl; 329 330 vtk << "DATASET UNSTRUCTURED_GRID\n"; 331 vtk << "POINTS " << bd.Nbubbles << " double" << std::endl; 332 for (size_t i = 0; i < bd.Nbubbles; ++i) 333 vtk << U[i][2] << " " << U[i][4] << " " << U[i][6] << std::endl; 334 vtk << std::endl; 335 336 vtk << "POINT_DATA " << bd.Nbubbles << std::endl; 337 vtk << "SCALARS radius double\n"; 338 vtk << "LOOKUP_TABLE default" << std::endl; 339 for (size_t i = 0; i < bd.Nbubbles; ++i) 340 vtk << U[i][0] << std::endl; 341 vtk << std::endl; 342 343 // vtk << "POINT_DATA " << bd.Nbubbles << std::endl; 344 vtk << "SCALARS velocity double\n"; 345 vtk << "LOOKUP_TABLE default" << std::endl; 346 for (size_t i = 0; i < bd.Nbubbles; ++i) 347 vtk << U[i][1] << std::endl; 348 vtk << std::endl; 349 350 vtk << "VECTORS transvel double\n"; 351 vtk << "LOOKUP_TABLE default" << std::endl; 352 for (size_t i = 0; i < bd.Nbubbles; ++i) 353 vtk << U[i][3] << " " << U[i][5] << " " << U[i][7] << std::endl; 354 vtk << std::endl; 355 } 356 } 357 }; 358 359 #endif /* BUBBLECLUSTERPOS_D_H_GHBNU569 */