bubble-dynamics

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

KellerMiksis.h (2865B)

      1 /* File:   KellerMiksis.h */
      2 /* Date:   Thu Feb 18 18:09:39 2016 */
      3 /* Author: Fabian Wermelinger */
      4 /* Tag:    Bubble dynamics based on Keller-Miksis (1980) */
      5 /* Copyright 2016 ETH Zurich. All Rights Reserved. */
      6 #ifndef KELLERMIKSIS_H_BMJWDX72
      7 #define KELLERMIKSIS_H_BMJWDX72
      8 
      9 #include "BubbleBase.h"
     10 
     11 #include <ODETB/GnuplotDump.h>
     12 #include <ODETB/TimeStepper/KernelBase.h>
     13 
     14 #include <cassert>
     15 #include <cmath>
     16 
     17 template <typename Tinput, typename Trhs=Tinput>
     18 class KellerMiksis : public KernelBase<Tinput,Trhs>
     19 {
     20     const size_t _N;
     21     GnuplotDump* m_dumper;
     22 
     23 public:
     24     KellerMiksis(const size_t N=0, const int format=GnuplotDump::BIN) : _N(N) { m_dumper = new GnuplotDump("out", format); }
     25     virtual ~KellerMiksis() { delete m_dumper;}
     26 
     27     static void IC(Tinput& U, BubbleData const& data)
     28     {
     29         assert(Tinput::DataType::SIZE == 2);
     30         assert(data.Nbubbles == 1);
     31         typename Tinput::DataType& IC = U[0];
     32         IC[0] = data.R0[0];
     33         IC[1] = data.Rdot0[0];
     34     }
     35 
     36     virtual void compute(const Tinput& U, Trhs& rhs, Real const t, void const* const data=nullptr)
     37     {
     38         const BubbleData& bd = *(BubbleData const* const)data;
     39         const typename Tinput::DataType& u = U[0];
     40         typename Trhs::DataType& r = rhs[0];
     41 
     42         assert(Tinput::DataType::SIZE == 2);
     43         assert(u[0] > 0.0);
     44         assert(u[1] < bd.cl);
     45         assert(data != nullptr);
     46         assert(bd.Nbubbles == 1);
     47 
     48         r[0] = u[1];
     49 
     50         const Real clInv = 1.0/bd.cl;
     51         const Real Rdc = u[1]*clInv;
     52         const Real pB = (bd.p0 - bd.pv + static_cast<Real>(2.0)*bd.sigma/bd.R0[0])*std::pow(bd.R0[0]/u[0], 3.0*bd.gamma); // bubble pressure
     53 
     54         r[1]  = (static_cast<Real>(0.5)*Rdc - static_cast<Real>(1.5))*u[1]*u[1];
     55         r[1] += bd.rInv*(static_cast<Real>(1.0) + static_cast<Real>(1.0 - 3.0*bd.gamma)*Rdc)*pB;
     56         r[1] -= bd.rInv*static_cast<Real>(2.0)*bd.sigma/u[0];
     57         r[1] -= static_cast<Real>(4.0)*bd.nuL*u[1]/u[0];
     58         r[1] -= bd.rInv*(static_cast<Real>(1.0) + Rdc)*(bd.p0 - bd.pv + bd.pext(t + u[0]*clInv));
     59         // r[1] -= bd.rInv*(static_cast<Real>(1.0) + Rdc)*(bd.p0 - bd.pv + bd.pext(t));
     60         r[1] -= bd.rInv*u[0]*clInv*bd.dpext(t + u[0]*clInv);
     61         // r[1] -= bd.rInv*u[0]*clInv*bd.dpext(t);
     62         r[1] /= (static_cast<Real>(1.0 - Rdc)*u[0] + static_cast<Real>(4.0)*bd.nuL*clInv);
     63 
     64         assert(!isnan(r[0]));
     65         assert(!isnan(r[1]));
     66     }
     67 
     68     virtual void write(const size_t step, const Real t, const Real dt, const Tinput& U, const void * const data=nullptr)
     69     {
     70         const BubbleData& bd = *(BubbleData const* const)data;
     71         const Real odat[3] = {U[0][0], U[0][1], bd.pext(t)};
     72         m_dumper->write(step, t, dt, &odat[0], 3);
     73         // m_dumper->write(step, t, dt, U.data(), Tinput::DataType::SIZE*U.size());
     74     }
     75 };
     76 
     77 #endif /* KELLERMIKSIS_H_BMJWDX72 */