bubble-dynamics

Spherical bubble dynamics simulator
git clone https://git.0xfab.ch/bubble-dynamics.git
Log | Files | Refs | README | LICENSE

KMCluster_TY.h (7383B)

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