commit f739cc31b8e32ca9ab1e080ac01eb804c583dac0
Author: Fabian Wermelinger <fabianw@mavt.ethz.ch>
Date: Thu, 22 Sep 2016 18:41:05 +0200
initial
Diffstat:
15 files changed, 961 insertions(+), 0 deletions(-)
diff --git a/.gitignore b/.gitignore
@@ -0,0 +1,5 @@
+*~
+*.bak
+*.swp
+*.swo
+bubbleDynamics
diff --git a/.gitmodules b/.gitmodules
@@ -0,0 +1,3 @@
+[submodule "ode_toolbox"]
+ path = ode_toolbox
+ url = scratch-petros.ethz.ch:/export/home/gitroot/ode_toolbox.git
diff --git a/Makefile b/Makefile
@@ -0,0 +1,65 @@
+CC = g++
+
+debug ?= 0
+omp ?= 0
+prec ?= double
+sse ?= 0
+align ?= 16
+sundials ?= 0
+eigen ?= 0
+sos = sg
+
+LIBS =
+CXXFLAGS = -std=c++11 -g -D_ALIGN_=$(align)
+CXXFLAGS += $(extra)
+
+ifeq "$(debug)" "0"
+ CXXFLAGS += -O3 -DNDEBUG
+ CXXFLAGS += -fno-expensive-optimizations -falign-functions=16
+else
+ CXXFLAGS += -O0
+endif
+ifeq "$(omp)" "1"
+ CXXFLAGS += -fopenmp -D_OMP_
+endif
+ifeq "$(prec)" "float"
+ CXXFLAGS += -D_FLOAT_PRECISION_
+ ifeq "$(sundials)" "1"
+ CXXFLAGS += -DSUNDIALS_SINGLE_PRECISION
+ endif
+else
+ CXXFLAGS += -D_DOUBLE_PRECISION_
+ ifeq "$(sundials)" "1"
+ CXXFLAGS += -DSUNDIALS_DOUBLE_PRECISION
+ endif
+endif
+ifeq "$(sse)" "1"
+ CXXFLAGS += -msse -msse2
+endif
+# ifeq "$(sos)" "prescribed" # default is sg (stiffened gas)
+# CXXFLAGS += -D_PRESCRIBED_SOS_
+# endif
+
+CXXFLAGS += -Iode_toolbox -Iode_toolbox/TimeStepper -Isrc
+ifeq "$(sundials)" "1"
+ CXXFLAGS += -D_USE_SUNDIALS_
+ CXXFLAGS += -IthirdParty/sundials/build/include
+ LIBS += -LthirdParty/sundials/build/lib -lsundials_cvode -lsundials_nvecserial
+endif
+ifeq "$(eigen)" "1"
+ CXXFLAGS += -D_USE_EIGEN_
+ CXXFLAGS += -I/usr/local/opt/eigen/include/eigen3
+endif
+
+VPATH := src
+SRC = src/bubbleDynamics.cpp
+HDR = $(wildcard *.h)
+
+.DEFAULT_GOAL := bubbleDynamics
+.PHONY := clean
+
+bubbleDynamics: $(SRC) $(HDR)
+ $(CC) $(CXXFLAGS) $(SRC) -o bubbleDynamics $(LIBS)
+
+clean:
+ rm -f *~ bubbleDynamics
diff --git a/examples/bremond2006/bremond.conf b/examples/bremond2006/bremond.conf
@@ -0,0 +1,32 @@
+# File : bremond.conf
+# Date : Thu Sep 22 13:15:27 2016
+# Author : Fabian Wermelinger
+# Description: Bubble cluster by Bremond et al. 2006
+# Copyright 2016 ETH Zurich. All Rights Reserved.
+
+# case selection
+dynamics kmcluster
+binary true
+
+# time stepper settings
+ts rkv56
+ts_tend 3.0e-5
+ts_dt 1.0e-9
+ts_writeGranularity 1000
+ts_reportGranularity 5000
+ts_maxLocalError 1.0e-7
+ts_rTol 1.0e-5
+ts_aTol 1.0e-5
+
+# bubble settings
+N 37
+bubbleFile bubblesBremont.dat
+rho 1000
+nu 0.0
+sigma 0.0
+gamma 1.4
+pv 0.0
+p0 1.0e5
+c 1481.0
+pamp 10.0
+pext coshbgl
diff --git a/examples/bremond2006/bubblesBremont.dat b/examples/bremond2006/bubblesBremont.dat
@@ -0,0 +1,38 @@
+37
+0.000000e+00 0.000000e+00 0.000000e+00 1.000000e-06 0.000000e+00
+2.000000e-04 0.000000e+00 0.000000e+00 1.000000e-06 0.000000e+00
+4.000000e-04 0.000000e+00 0.000000e+00 1.000000e-06 0.000000e+00
+6.000000e-04 0.000000e+00 0.000000e+00 1.000000e-06 0.000000e+00
+3.000000e-04 1.732051e-04 0.000000e+00 1.000000e-06 0.000000e+00
+5.000000e-04 1.732051e-04 0.000000e+00 1.000000e-06 0.000000e+00
+4.000000e-04 3.464102e-04 0.000000e+00 1.000000e-06 0.000000e+00
+1.000000e-04 1.732051e-04 0.000000e+00 1.000000e-06 0.000000e+00
+2.000000e-04 3.464102e-04 0.000000e+00 1.000000e-06 0.000000e+00
+3.000000e-04 5.196152e-04 0.000000e+00 1.000000e-06 0.000000e+00
+5.421011e-20 3.464102e-04 0.000000e+00 1.000000e-06 0.000000e+00
+1.000000e-04 5.196152e-04 0.000000e+00 1.000000e-06 0.000000e+00
+-1.000000e-04 5.196152e-04 0.000000e+00 1.000000e-06 0.000000e+00
+-1.000000e-04 1.732051e-04 0.000000e+00 1.000000e-06 0.000000e+00
+-2.000000e-04 3.464102e-04 0.000000e+00 1.000000e-06 0.000000e+00
+-3.000000e-04 5.196152e-04 0.000000e+00 1.000000e-06 0.000000e+00
+-3.000000e-04 1.732051e-04 0.000000e+00 1.000000e-06 0.000000e+00
+-4.000000e-04 3.464102e-04 0.000000e+00 1.000000e-06 0.000000e+00
+-5.000000e-04 1.732051e-04 0.000000e+00 1.000000e-06 0.000000e+00
+-2.000000e-04 2.449294e-20 0.000000e+00 1.000000e-06 0.000000e+00
+-4.000000e-04 4.898587e-20 0.000000e+00 1.000000e-06 0.000000e+00
+-6.000000e-04 7.347881e-20 0.000000e+00 1.000000e-06 0.000000e+00
+-3.000000e-04 -1.732051e-04 0.000000e+00 1.000000e-06 0.000000e+00
+-5.000000e-04 -1.732051e-04 0.000000e+00 1.000000e-06 0.000000e+00
+-4.000000e-04 -3.464102e-04 0.000000e+00 1.000000e-06 0.000000e+00
+-1.000000e-04 -1.732051e-04 0.000000e+00 1.000000e-06 0.000000e+00
+-2.000000e-04 -3.464102e-04 0.000000e+00 1.000000e-06 0.000000e+00
+-3.000000e-04 -5.196152e-04 0.000000e+00 1.000000e-06 0.000000e+00
+-1.897354e-19 -3.464102e-04 0.000000e+00 1.000000e-06 0.000000e+00
+-1.000000e-04 -5.196152e-04 0.000000e+00 1.000000e-06 0.000000e+00
+1.000000e-04 -5.196152e-04 0.000000e+00 1.000000e-06 0.000000e+00
+1.000000e-04 -1.732051e-04 0.000000e+00 1.000000e-06 0.000000e+00
+2.000000e-04 -3.464102e-04 0.000000e+00 1.000000e-06 0.000000e+00
+3.000000e-04 -5.196152e-04 0.000000e+00 1.000000e-06 0.000000e+00
+3.000000e-04 -1.732051e-04 0.000000e+00 1.000000e-06 0.000000e+00
+4.000000e-04 -3.464102e-04 0.000000e+00 1.000000e-06 0.000000e+00
+5.000000e-04 -1.732051e-04 0.000000e+00 1.000000e-06 0.000000e+00
diff --git a/examples/keller-miksis/km.conf b/examples/keller-miksis/km.conf
@@ -0,0 +1,30 @@
+# File : km.conf
+# Date : Thu Sep 22 13:15:27 2016
+# Author : Fabian Wermelinger
+# Description: Keller-Miksis example
+# Copyright 2016 ETH Zurich. All Rights Reserved.
+
+# case selection
+dynamics km # Keller-Miksis
+
+# time stepper settings
+ts rkv56
+ts_tend 3.0e-5
+ts_dt 1.0e-12
+ts_writeGranularity 1
+ts_reportGranularity 5000
+ts_maxLocalError 1.0e-8
+ts_rTol 1.0e-6
+ts_aTol 1.0e-6
+
+# bubble settings
+R0 4.0e-6
+rho 1000
+nu 0.0
+sigma 0.0
+gamma 1.4
+pv 0.0
+p0 1.0e5
+c 1481.0
+pamp 2.4
+pext coshbgl
diff --git a/ode_toolbox b/ode_toolbox
@@ -0,0 +1 @@
+Subproject commit 1bc8a2d5cbaf162ce4f10aa90342b7bbc3008f5e
diff --git a/src/BubbleBase.h b/src/BubbleBase.h
@@ -0,0 +1,150 @@
+/* File: BubbleBase.h */
+/* Date: Fri Feb 12 12:39:08 2016 */
+/* Author: Fabian Wermelinger */
+/* Tag: Bubble base class */
+/* Copyright 2016 ETH Zurich. All Rights Reserved. */
+#ifndef BUBBLEBASE_H_ZVSZD3E8
+#define BUBBLEBASE_H_ZVSZD3E8
+
+#include <vector>
+#include <string>
+#include <cmath>
+#include <cstdio>
+#include <functional>
+#include "ArgumentParser.h"
+
+struct Bubble
+{
+ double pos[3];
+ double distance(const Bubble& b) const
+ {
+ const double dx = b.pos[0] - pos[0];
+ const double dy = b.pos[1] - pos[1];
+ const double dz = b.pos[2] - pos[2];
+ return std::sqrt(dx*dx + dy*dy + dz*dz);
+ }
+};
+
+struct BubbleData
+{
+ BubbleData(ArgumentParser& p)
+ {
+ Nbubbles = p("-N").asInt(1);
+ R0.reserve(Nbubbles);
+ Rdot0.reserve(Nbubbles);
+ coords.reserve(Nbubbles);
+ Dij.reserve(Nbubbles*Nbubbles);
+
+ R0[0] = p("-R0").asDouble(4.0e-6); // m
+ Rdot0[0] = p("-Rdot0").asDouble(0.0); // m/s
+ coords[0].pos[0] = coords[0].pos[1] = coords[0].pos[2] = 0.0;
+
+ // gas
+ gamma = p("-gamma").asDouble(1.4);
+
+ // liquid
+ rhoL = p("-rho").asDouble(1000.0); // kg/m^3
+ nuL = p("-nu").asDouble(1.0e-6); // m^2/s
+ pv = p("-pv").asDouble(2400.0); // Pa
+ p0 = p("-p0").asDouble(1.0e5); // Pa
+ sigma = p("-sigma").asDouble(7.3e-2); // kg/s^2
+ h = p("-h").asDouble(4.514673e-7);
+ cl = p("-c").asDouble(1481.0); // m/s
+
+ // external forcing
+ pamp = p("-pamp").asDouble(1.4); // external pressure amplification
+ omega = p("-omega").asDouble(2.0*M_PI*26500.0); // s^-1
+ p_pre = p("-p_pre").asDouble(1.0e5);
+ p_post= p("-p_post").asDouble(1.0e5);
+ t_12 = p("-t_12").asDouble(1.0e-6);
+ dt_smooth = p("-dt_smooth").asDouble(1.0e-9);
+
+ // drivers
+ if (p("-pext").asString("const") == "const")
+ pext = [this] (const Real t) { return (this->pamp * this->p0); };
+ if (p("-pext").asString("const") == "sin")
+ pext = [this] (const Real t) { return (-this->p0*this->pamp*std::sin(this->omega * t)); };
+ if (p("-pext").asString("const") == "coshbgl")
+ pext = [this] (const Real t) { return (this->p0 * (1.0 - this->pamp*std::cos(this->omega * (t - t_star)))); };
+ if (p("-pext").asString("const") == "cosramp")
+ pext = [this] (const Real t)
+ {
+ const double phi = std::min(1.0,std::max(0.0,(t-this->t_12)/this->dt_smooth+0.5));
+ return this->p_post + 0.5*(this->p_pre-this->p_post)*(1.0+std::cos(M_PI*phi));
+ };
+
+ // pressure derivative
+ dpext = [this] (const Real t) { return 0.0; };
+
+ // other
+ filename = p("-bubbleFile").asString("none");
+ bVTK = p("-vtk").asBool(false);
+
+ // helpers
+ assert(rhoL > 0);
+ rInv = 1.0/rhoL;
+ t_star = 1.0/omega * std::acos(1.0/pamp);
+ }
+
+ // Number of bubbles and their initial state
+ size_t Nbubbles;
+ std::vector<Real> R0;
+ std::vector<Real> Rdot0;
+ std::vector<Bubble> coords;
+
+ // Bubble properties
+ Real gamma; // ratio of specific heats for gas
+
+ // Liquid properties
+ Real rhoL; // density
+ Real nuL; // viscosity
+ Real pv; // vaporization pressure
+ Real p0; // far field pressure
+ Real sigma; // surface tension
+ Real h; // van der Waals hard-core radius
+ Real cl; // speed of sound
+
+ // external forcing
+ Real pamp; // forcing pressure amplitude
+ Real omega; // external frequency
+ Real p_pre, p_post;
+ Real t_12, dt_smooth;
+
+ // driver functions
+ std::function<Real(const Real)> pext;
+ std::function<Real(const Real)> dpext;
+
+ // helpers
+ std::string filename;
+ std::vector<Real> Dij;
+ Real rInv;
+ Real t_star;
+ bool bVTK;
+
+ template <typename T>
+ inline Real pBubble(const T& U, const size_t i=0) const
+ {
+ return (p0 - pv + 2.0*sigma/R0[i])*std::pow(R0[i]/U[i][0], 3.0*gamma) - 2.0*sigma/U[i][0] - 4.0*nuL*rhoL/U[i][0]*U[i][1];
+ }
+
+ void print() const
+ {
+ printf("Bubble Data:\n");
+ printf("\tNbubbles = %d\n", Nbubbles);
+ printf("\tBubbles:\n");
+ for (size_t i = 0; i < Nbubbles; ++i)
+ printf("\t\tBubble %d: R0 = %e, Rdot0 = %e\n", i, R0[i], Rdot0[i]);
+ printf("\tgamma = %e\n", gamma);
+ printf("\trho = %e\n", rhoL);
+ printf("\tnu = %e\n", nuL);
+ printf("\tpv = %e\n", pv);
+ printf("\tp0 = %e\n", p0);
+ printf("\tsigma = %e\n", sigma);
+ printf("\th = %e\n", h);
+ printf("\tcl = %e\n", cl);
+ printf("\tpamp = %e\n", pamp);
+ printf("\tomega = %e\n", omega);
+ }
+};
+
+#endif /* BUBBLEBASE_H_ZVSZD3E8 */
diff --git a/src/bubbleDynamics.cpp b/src/bubbleDynamics.cpp
@@ -0,0 +1,150 @@
+/* File: bubbleDynamics.cpp */
+/* Date: Tue Jan 12 22:04:58 2016 */
+/* Author: Fabian Wermelinger */
+/* Tag: Solve 1D bubble dynamics */
+/* Copyright 2016 ETH Zurich. All Rights Reserved. */
+#include <cassert>
+#include <cstdio>
+#include <iostream>
+#include <iomanip>
+#include <fstream>
+#include <cmath>
+#include <ctime>
+
+#include "common.h"
+#include "kernels.h"
+#include "GnuplotDump.h"
+#include "TimeStepper.h"
+
+using namespace std;
+
+
+void say_hello(const int argc, const char * argv[])
+{
+ time_t rawtime;
+ struct tm* timeinfo;
+ time(&rawtime);
+ timeinfo = localtime(&rawtime);
+ char buf[256];
+ strftime(buf, 256, "%A, %h %d %Y, %r", timeinfo);
+ cout << "Hi! It is now " << buf << endl;
+ cout << "Here is what you gave me: ";
+ for (int i = 0; i < argc; ++i)
+ cout << argv[i] << " ";
+ cout << endl;
+}
+
+template <typename Tinput, typename Trhs>
+KernelBase<Tinput,Trhs>* kernelFactory(ArgumentParser& parser, Tinput& U, BubbleData& simConfig)
+{
+ int yeaBinary;
+ if (parser("-binary").asBool(false))
+ yeaBinary = GnuplotDump::BIN;
+ else
+ yeaBinary = GnuplotDump::ASCII;
+
+ if (parser("-dynamics").asString("rp") == "rp")
+ {
+ RayleighPlesset<Tinput,Trhs>::IC(U, simConfig);
+ return new RayleighPlesset<Tinput,Trhs>(U.size(), yeaBinary);
+ }
+ else if (parser("-dynamics").asString("rp") == "hbgl")
+ {
+ RayleighPlesset_HBGL<Tinput,Trhs>::IC(U, simConfig);
+ return new RayleighPlesset_HBGL<Tinput,Trhs>(U.size(), yeaBinary);
+ }
+ else if (parser("-dynamics").asString("rp") == "km")
+ {
+ KellerMiksis<Tinput,Trhs>::IC(U, simConfig);
+ return new KellerMiksis<Tinput,Trhs>(U.size(), yeaBinary);
+ }
+#ifdef _USE_EIGEN_
+ else if (parser("-dynamics").asString("rp") == "kmcluster")
+ {
+ KMCluster_TY<Tinput,Trhs>::IC(U, simConfig);
+ return new KMCluster_TY<Tinput,Trhs>(U.size(), yeaBinary);
+ }
+#endif
+ else
+ return nullptr;
+}
+
+template <typename Tinput, typename Trhs>
+StepperBase<Tinput,Trhs>* stepperFactory(ArgumentParser& p, Tinput& U, StepperSettings& S, BubbleData& D)
+{
+ KernelBase<Tinput,Trhs> * const kern = kernelFactory<Tinput,Trhs>(p, U, D);
+ assert(kern != nullptr);
+ return S.template stepperFactory<Tinput,Trhs>(p, U, kern);
+}
+
+
+int main(int argc, const char** argv)
+{
+ say_hello(argc, argv);
+
+ // Parameter
+ ArgumentParser parser(argc, argv);
+ if (parser.exist("-conf")) parser.readFile(parser("-conf").asString());
+ parser.print_args();
+
+ // set user specific data
+ StepperSettings stepperData(parser);
+ BubbleData simConfig(parser);
+
+ // set up solution vector
+ using vec_t = StateVector<State2>;
+ vec_t U(simConfig.Nbubbles);
+
+ // Select stepper and kernel
+ StepperBase<vec_t,vec_t> * const stepper = stepperFactory<vec_t,vec_t>(parser, U, stepperData, simConfig);
+ assert(stepper != nullptr);
+
+ simConfig.print();
+
+ size_t& step = stepperData.step;
+ double& t = stepperData.t;
+ double& dt = stepperData.dt;
+ const double& tFinal = stepperData.tFinal;
+ double& tDump = stepperData.tDump;
+ if (stepperData.dtDump > 0.0) tDump = t + stepperData.dtDump;
+
+ stepper->write(step, t, dt, U, &simConfig);
+
+ while (t < tFinal)
+ {
+ dt = (tFinal-t) < dt ? (tFinal-t) : dt;
+ if (stepperData.dtDump > 0.0) dt = (tDump-t) < dt ? (tDump-t) : dt;
+
+ stepper->step(&simConfig);
+
+ if (stepperData.bFixedStep && stepperData.dtDump > 0.0)
+ {
+ if (t == tDump)
+ {
+ stepper->write(step, t, dt, U, &simConfig);
+ tDump += stepperData.dtDump;
+ }
+ }
+ else
+ if (stepperData.bFixedStep && step % stepperData.writeGranularity == 0)
+ stepper->write(step, t, dt, U, &simConfig);
+ }
+
+ if (stepperData.dtDump > 0.0)
+ {
+ if (t != tDump)
+ stepper->write(step, t, dt, U, &simConfig);
+ }
+ else
+ {
+ if (step % stepperData.writeGranularity != 0)
+ stepper->write(step, t, dt, U, &simConfig);
+ }
+
+ printf("\nFinal Time = %e;\nFinal Step = %d\n", stepperData.t, stepperData.step);
+
+ stepper->dispose();
+ delete stepper;
+
+ return 0;
+}
diff --git a/src/kernels.h b/src/kernels.h
@@ -0,0 +1,16 @@
+/* File: kernels.h */
+/* Date: Tue Jan 12 22:57:38 2016 */
+/* Author: Fabian Wermelinger */
+/* Tag: Right-hand side kernels */
+/* Copyright 2016 ETH Zurich. All Rights Reserved. */
+#ifndef KERNELS_H_KY4W86QV
+#define KERNELS_H_KY4W86QV
+
+#include "kernels/RayleighPlesset.h"
+#include "kernels/RayleighPlesset_HBGL.h"
+#include "kernels/KellerMiksis.h"
+#ifdef _USE_EIGEN_
+#include "kernels/KMCluster_TY.h"
+#endif
+
+#endif /* KERNELS_H_KY4W86QV */
diff --git a/src/kernels/Dahlquist.h b/src/kernels/Dahlquist.h
@@ -0,0 +1,56 @@
+/* File: Dahlquist.h */
+/* Date: Fri Feb 5 19:31:02 2016 */
+/* Author: Fabian Wermelinger */
+/* Tag: Dahlquist equation (Test kernel) */
+/* Copyright 2016 ETH Zurich. All Rights Reserved. */
+#ifndef DAHLQUIST_H_SQ60WZPY
+#define DAHLQUIST_H_SQ60WZPY
+
+#include <cassert>
+#include "common.h"
+#include "KernelBase.h"
+#include "GnuplotDump.h"
+
+template <typename Tinput, typename Trhs>
+class Dahlquist : public KernelBase<Tinput,Trhs>
+{
+ GnuplotDump m_gp;
+
+protected:
+ const Real m_a;
+
+public:
+ Dahlquist(const Real a=1.0) : m_gp("out", GnuplotDump::ASCII), m_a(a) {}
+ ~Dahlquist() {}
+
+ void compute(const Tinput& U, Trhs& rhs, const Real t, void const* const data=nullptr)
+ {
+ rhs = m_a*U;
+ }
+
+ void write(const size_t step, const Real t, const Real dt, const Tinput& U, const void * const data=nullptr)
+ {
+ m_gp.write(step, t, dt, U.data(), U.size()*Tinput::DataType::SIZE);
+ }
+};
+
+template <typename Tinput, typename Trhs>
+class DahlquistLSRK3 : public Dahlquist<Tinput,Trhs>
+{
+ using Dahlquist<Tinput, Trhs>::m_a;
+
+public:
+ DahlquistLSRK3(const Real a=1.0) : Dahlquist<Tinput,Trhs>(a) {}
+ ~DahlquistLSRK3() {}
+
+ void compute(const Tinput& U, Trhs& rhs, const Real t, void const* const data=nullptr)
+ {
+ assert(data != nullptr);
+ const Real * const LSRKData = static_cast<const Real* const>(data);
+ const Real dt = LSRKData[0];
+ const Real A = LSRKData[1];
+ rhs = A*rhs + dt*m_a*U;
+ }
+};
+
+#endif /* DAHLQUIST_H_SQ60WZPY */
diff --git a/src/kernels/KMCluster_TY.h b/src/kernels/KMCluster_TY.h
@@ -0,0 +1,208 @@
+/* File: KMCluster_TY.h */
+/* Date: Mon Feb 22 14:48:12 2016 */
+/* Author: Fabian Wermelinger */
+/* Tag: Bubble cluster based on Keller-Miksis following Tiwari (2015) */
+/* Copyright 2016 ETH Zurich. All Rights Reserved. */
+#ifndef BUBBLECLUSTER_TY_H_ECCX59FB
+#define BUBBLECLUSTER_TY_H_ECCX59FB
+
+#include <cassert>
+#include <cstdlib>
+#include <fstream>
+#include <cmath>
+#include <vector>
+
+#include <Eigen/Core>
+#include <Eigen/Dense>
+
+#include "KernelBase.h"
+#include "BubbleBase.h"
+#include "GnuplotDump.h"
+
+
+template <typename Tinput, typename Trhs>
+class KMCluster_TY : public KernelBase<Tinput,Trhs>
+{
+ const size_t _N;
+ GnuplotDump* m_dumper;
+
+public:
+ KMCluster_TY(const size_t N=0, const int format=GnuplotDump::BIN) : _N(N) { m_dumper = new GnuplotDump("out", format); }
+ virtual ~KMCluster_TY() { delete m_dumper;}
+
+ static void IC(Tinput& U, BubbleData& data)
+ {
+ assert(Tinput::DataType::SIZE == 2);
+ assert(data.Nbubbles == U.size());
+ const size_t N = data.Nbubbles;
+
+ // read IC from file
+ std::ifstream infile(data.filename);
+ size_t k;
+ infile >> k;
+ if (N != k)
+ {
+ std::cerr << "KMCluster_TY: Allocated " << N << " bubbles."<< std::endl;
+ std::cerr << " File " << data.filename << " expects " << k << " bubbles."<< std::endl;
+ abort();
+ }
+
+ for (size_t i = 0; i < N; ++i)
+ {
+ if (infile.good())
+ {
+ Bubble& b = data.coords[i];
+ infile >> b.pos[0] >> b.pos[1] >> b.pos[2];
+ infile >> data.R0[i];
+ infile >> data.Rdot0[i];
+
+ State2& IC = U[i];
+ IC[0] = data.R0[i];
+ IC[1] = data.Rdot0[i];
+ }
+ else
+ {
+ std::cerr << "KMCluster_TY: Input file has incorrect format" << std::endl;
+ abort();
+ }
+ }
+
+ for (int i = 0; i < N; ++i)
+ {
+ const Bubble& bi = data.coords[i];
+ for (int j = 0; j < N; ++j)
+ {
+ const Bubble& bj = data.coords[j];
+ data.Dij[j + i*N] = bi.distance(bj);
+ }
+ }
+
+#ifndef NDEBUG
+ for (size_t i = 0; i < N; ++i)
+ assert(0 == data.Dij[i + i*N]);
+#endif /* NDEBUG */
+ }
+
+ virtual void compute(const Tinput& U, Trhs& rhs, Real const t, void const* const data=nullptr)
+ {
+ const BubbleData& bd = *(BubbleData const* const)data;
+
+ assert(_N == U.size());
+ Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic> A(_N,_N);
+ Eigen::Matrix<Real, Eigen::Dynamic, 1> b(_N);
+
+ const Real clInv = 1.0/bd.cl;
+
+ for (size_t i = 0; i < _N; ++i)
+ {
+ const State2& Ui = U[i];
+ // fill row i of A
+ Real bnbr = 0;
+ for (size_t j = 0; j < _N; ++j)
+ {
+ const State2& Uj = U[j];
+ if (i == j)
+ A(i,j) = (static_cast<Real>(1) - Ui[1]*clInv)*Ui[0] + static_cast<Real>(4.0)*bd.nuL*clInv;
+ else
+ {
+ const Real UDij = Uj[0]/bd.Dij[j + i*_N];
+ A(i,j) = Uj[0]*UDij;
+ bnbr += Uj[1]*Uj[1]*UDij;
+ }
+ }
+
+ // compute bi
+ const Real pstat = bd.p0 - bd.pv;
+ const Real Rdc = Ui[1]*clInv;
+ 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);
+
+ b(i) = (static_cast<Real>(0.5)*Rdc - static_cast<Real>(1.5))*Ui[1]*Ui[1];
+ b(i) -= bd.rInv*(static_cast<Real>(1.0) + Rdc)*(pstat + bd.pext(t + Ui[0]*clInv));
+ // b(i) -= bd.rInv*(static_cast<Real>(1.0) + Rdc)*(pstat + bd.pext(t));
+ 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];
+ b(i) += bd.rInv*static_cast<Real>(1.0 - 3.0*bd.gamma)*Rdc*pbs;
+ b(i) -= bd.rInv*Rdc*bd.dpext(t + Ui[0]*clInv);
+ // b(i) -= bd.rInv*Rdc*bd.dpext(t);
+ b(i) -= static_cast<Real>(2.0)*bnbr;
+
+ assert(!isnan(b(i)));
+ }
+
+ /* std::cout << A << std::endl << std::endl; */
+ Eigen::Matrix<Real, Eigen::Dynamic, 1> Rddot = A.fullPivHouseholderQr().solve(b);
+
+ for (size_t i = 0; i < _N; ++i)
+ {
+ const State2& Ui = U[i];
+ State2& ri = rhs[i];
+ ri[0] = Ui[1];
+ ri[1] = Rddot(i);
+ }
+ }
+
+ virtual void write(const size_t step, const Real t, const Real dt, const Tinput& U, const void * const data=nullptr)
+ {
+ static size_t callCount = 0;
+
+ const BubbleData& bd = *(BubbleData const* const)data;
+
+ // double pb[3];
+ // const int id[3] = {0, 7, 1};
+ // for (int i = 0; i < 3; ++i)
+ // pb[i] = bd.pBubble(U, id[i]);
+
+ // // class 1, 2, 3, 4, 5, 6
+ // // const Real odat[6] = {U[0][0], U[1][0], U[4][0], U[2][0], U[5][0], U[3][0]};
+ // const Real odat[7] = {U[0][0], U[7][0], U[1][0], pb[0], pb[1], pb[2], bd.pext(t)};
+
+ const size_t N = 102;
+ double out[N];
+ double Vg = 0.0;
+ for (int i = 0; i < 50; ++i)
+ {
+ const double r = U[i][0];
+ out[i] = r;
+ out[i+50] = bd.pBubble(U, i);
+ Vg += 4.0/3.0*r*r*r*M_PI;
+ }
+ out[100] = bd.pext(t);
+ out[101] = Vg;
+
+ m_dumper->write(step, t, dt, &out[0], N);
+
+ if (bd.bVTK)
+ {
+ // dump bubble vtk
+ std::stringstream streamer;
+ streamer<< "bubbles-vtk_" << setfill('0') << setw(8) << callCount++ << ".vtk";
+ std::ofstream vtk(streamer.str().c_str());
+ vtk.setf(std::ios::scientific, std::ios::floatfield);
+
+ vtk << "# vtk DataFile Version 2.0\n";
+ vtk << "Unstructured grid legacy vtk file with point scalar data\n";
+ vtk << "ASCII\n" << std::endl;
+
+ vtk << "DATASET UNSTRUCTURED_GRID\n";
+ vtk << "POINTS " << bd.Nbubbles << " double" << std::endl;
+ for (int i = 0; i < bd.Nbubbles; ++i)
+ vtk << bd.coords[i].pos[0] << " " << bd.coords[i].pos[1] << " " << bd.coords[i].pos[2] << std::endl;
+ vtk << std::endl;
+
+ vtk << "POINT_DATA " << bd.Nbubbles << std::endl;
+ vtk << "SCALARS radius double\n";
+ vtk << "LOOKUP_TABLE default" << std::endl;
+ for (int i = 0; i < bd.Nbubbles; ++i)
+ vtk << U[i][0] << std::endl;
+ vtk << std::endl;
+
+ // vtk << "POINT_DATA " << bd.Nbubbles << std::endl;
+ vtk << "SCALARS velocity double\n";
+ vtk << "LOOKUP_TABLE default" << std::endl;
+ for (int i = 0; i < bd.Nbubbles; ++i)
+ vtk << U[i][1] << std::endl;
+ vtk << std::endl;
+ }
+ }
+};
+
+#endif /* BUBBLECLUSTER_TY_H_ECCX59FB */
diff --git a/src/kernels/KellerMiksis.h b/src/kernels/KellerMiksis.h
@@ -0,0 +1,74 @@
+/* File: KellerMiksis.h */
+/* Date: Thu Feb 18 18:09:39 2016 */
+/* Author: Fabian Wermelinger */
+/* Tag: Bubble dynamics based on Keller-Miksis (1980) */
+/* Copyright 2016 ETH Zurich. All Rights Reserved. */
+#ifndef KELLERMIKSIS_H_BMJWDX72
+#define KELLERMIKSIS_H_BMJWDX72
+
+#include <cassert>
+#include <cmath>
+
+#include "KernelBase.h"
+#include "BubbleBase.h"
+#include "GnuplotDump.h"
+
+template <typename Tinput, typename Trhs>
+class KellerMiksis : public KernelBase<Tinput,Trhs>
+{
+ const size_t _N;
+ GnuplotDump* m_dumper;
+
+public:
+ KellerMiksis(const size_t N=0, const int format=GnuplotDump::BIN) : _N(N) { m_dumper = new GnuplotDump("out", format); }
+ virtual ~KellerMiksis() { delete m_dumper;}
+
+ static void IC(Tinput& U, BubbleData const& data)
+ {
+ assert(Tinput::DataType::SIZE == 2);
+ assert(data.Nbubbles == 1);
+ State2& IC = U[0];
+ IC[0] = data.R0[0];
+ IC[1] = data.Rdot0[0];
+ }
+
+ virtual void compute(const Tinput& U, Trhs& rhs, Real const t, void const* const data=nullptr)
+ {
+ const BubbleData& bd = *(BubbleData const* const)data;
+ const State2& u = U[0];
+ State2& r = rhs[0];
+
+ assert(Tinput::DataType::SIZE == 2);
+ assert(u[0] > 0.0);
+ assert(u[1] < bd.cl);
+ assert(data != nullptr);
+ assert(bd.Nbubbles == 1);
+
+ r[0] = u[1];
+
+ const Real pstat = bd.p0 - bd.pv;
+ const Real Rdc = u[1]/bd.cl;
+ const Real pbs = (pstat + static_cast<Real>(2.0)*bd.sigma/bd.R0[0])*std::pow(bd.R0[0]/u[0], 3.0*bd.gamma);
+ r[1] = (static_cast<Real>(0.5)*Rdc - static_cast<Real>(1.5))*u[1]*u[1];
+ r[1] -= bd.rInv*(static_cast<Real>(1.0) + Rdc)*(pstat + bd.pext(t + u[0]/bd.cl));
+ // r[1] -= bd.rInv*(static_cast<Real>(1.0) + Rdc)*(pstat + bd.pext(t));
+ r[1] += bd.rInv*(pbs - static_cast<Real>(2.0)*bd.sigma/u[0]) - static_cast<Real>(4.0)*bd.nuL*u[1]/u[0];
+ r[1] += bd.rInv*static_cast<Real>(1.0 - 3.0*bd.gamma)*Rdc*pbs;
+ r[1] -= bd.rInv*Rdc*bd.dpext(t + u[0]/bd.cl);
+ // r[1] -= bd.rInv*Rdc*bd.dpext(t);
+ r[1] /= (static_cast<Real>(1.0 - Rdc)*u[0] + static_cast<Real>(4.0)*bd.nuL/bd.cl);
+
+ assert(!isnan(r[0]));
+ assert(!isnan(r[1]));
+ }
+
+ virtual void write(const size_t step, const Real t, const Real dt, const Tinput& U, const void * const data=nullptr)
+ {
+ const BubbleData& bd = *(BubbleData const* const)data;
+ const Real odat[3] = {U[0][0], U[0][1], bd.pext(t)};
+ m_dumper->write(step, t, dt, &odat[0], 3);
+ // m_dumper->write(step, t, dt, U.data(), Tinput::DataType::SIZE*U.size());
+ }
+};
+
+#endif /* KELLERMIKSIS_H_BMJWDX72 */
diff --git a/src/kernels/RayleighPlesset.h b/src/kernels/RayleighPlesset.h
@@ -0,0 +1,65 @@
+/* File: RayleighPlesset.h */
+/* Date: Tue Jan 12 23:00:47 2016 */
+/* Author: Fabian Wermelinger */
+/* Tag: Rayleigh-Plesset dynamics */
+/* Copyright 2016 ETH Zurich. All Rights Reserved. */
+#ifndef RAYLEIGHPLESSET_H_DLCI4XQ3
+#define RAYLEIGHPLESSET_H_DLCI4XQ3
+
+#include <cassert>
+#include <cmath>
+
+#include "KernelBase.h"
+#include "BubbleBase.h"
+#include "GnuplotDump.h"
+
+template <typename Tinput, typename Trhs>
+class RayleighPlesset : public KernelBase<Tinput,Trhs>
+{
+ const size_t _N;
+ GnuplotDump* m_dumper;
+
+public:
+ RayleighPlesset(const size_t N=0, const int format=GnuplotDump::BIN) : _N(N) { m_dumper = new GnuplotDump("out", format); }
+ virtual ~RayleighPlesset() { delete m_dumper;}
+
+ static void IC(Tinput& U, BubbleData const& data)
+ {
+ assert(Tinput::DataType::SIZE == 2);
+ assert(data.Nbubbles == 1);
+ State2& IC = U[0];
+ IC[0] = data.R0[0];
+ IC[1] = data.Rdot0[0];
+ }
+
+ virtual void compute(const Tinput& U, Trhs& rhs, Real const t, void const* const data=nullptr)
+ {
+ const BubbleData& bd = *(BubbleData const* const)data;
+ const State2& u = U[0];
+ State2& r = rhs[0];
+
+ assert(Tinput::DataType::SIZE == 2);
+ assert(u[0] > 0.0);
+ assert(data != nullptr);
+ assert(bd.Nbubbles == 1);
+
+ r[0] = u[1];
+
+ r[1] = bd.rInv*(bd.pv - bd.pext(t));
+ r[1] += bd.rInv*(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);
+ r[1] -= static_cast<Real>(1.5)*u[1]*u[1];
+ r[1] -= static_cast<Real>(4.0)*bd.nuL*u[1]/u[0];
+ r[1] -= static_cast<Real>(2.0)*bd.sigma*bd.rInv/u[0];
+ r[1] /= u[0];
+
+ assert(!isnan(r[0]));
+ assert(!isnan(r[1]));
+ }
+
+ virtual void write(const size_t step, const Real t, const Real dt, const Tinput& U, const void * const data=nullptr)
+ {
+ m_dumper->write(step, t, dt, U.data(), Tinput::DataType::SIZE*U.size());
+ }
+};
+
+#endif /* RAYLEIGHPLESSET_H_DLCI4XQ3 */
diff --git a/src/kernels/RayleighPlesset_HBGL.h b/src/kernels/RayleighPlesset_HBGL.h
@@ -0,0 +1,68 @@
+/* File: RayleighPlesset_HBGL.h */
+/* Date: Thu Feb 18 18:23:07 2016 */
+/* Author: Fabian Wermelinger */
+/* Tag: Rayleigh-Plesset dynamics (Hilgenfeldt, Brenner, Grossmann, Lohse, 1998) */
+/* Copyright 2016 ETH Zurich. All Rights Reserved. */
+#ifndef RAYLEIGHPLESSET_HBGL_H_EVWDWXSM
+#define RAYLEIGHPLESSET_HBGL_H_EVWDWXSM
+
+#include <cassert>
+#include <cmath>
+
+#include "KernelBase.h"
+#include "BubbleBase.h"
+#include "GnuplotDump.h"
+
+template <typename Tinput, typename Trhs>
+class RayleighPlesset_HBGL : public KernelBase<Tinput,Trhs>
+{
+ const size_t _N;
+ GnuplotDump* m_dumper;
+
+public:
+ RayleighPlesset_HBGL(const size_t N=0, const int format=GnuplotDump::BIN) : _N(N) { m_dumper = new GnuplotDump("out", format); }
+ virtual ~RayleighPlesset_HBGL() { delete m_dumper;}
+
+ static void IC(Tinput& U, BubbleData const& data)
+ {
+ assert(Tinput::DataType::SIZE == 2);
+ assert(data.Nbubbles == 1);
+ State2& IC = U[0];
+ IC[0] = data.R0[0];
+ IC[1] = data.Rdot0[0];
+ }
+
+ virtual void compute(const Tinput& U, Trhs& rhs, Real const t, void const* const data=nullptr)
+ {
+ const BubbleData& bd = *(BubbleData const* const)data;
+ const State2& u = U[0];
+ State2& r = rhs[0];
+
+ assert(Tinput::DataType::SIZE == 2);
+ assert(u[0] > 0.0);
+ assert(u[1] < bd.cl);
+ assert(data != nullptr);
+ assert(bd.Nbubbles == 1);
+
+ r[0] = u[1];
+
+ const Real dR03 = (bd.R0[0]*bd.R0[0]*bd.R0[0] - bd.h*bd.h*bd.h);
+ const Real dR3 = (u[0]*u[0]*u[0] - bd.h*bd.h*bd.h);
+ r[1] = bd.rInv*(bd.p0 + static_cast<Real>(2.0)*bd.sigma/bd.R0[0])*std::pow(dR03/dR3,bd.gamma)*(static_cast<Real>(1.0) - static_cast<Real>(3.0)*bd.gamma*u[0]*u[0]*u[0]*u[1]/(bd.cl*dR3));
+ r[1] -= bd.rInv*bd.pext(t);
+ r[1] -= static_cast<Real>(1.5)*u[1]*u[1];
+ r[1] -= static_cast<Real>(4.0)*bd.nuL*u[1]/u[0];
+ r[1] -= static_cast<Real>(2.0)*bd.sigma*bd.rInv/u[0];
+ r[1] /= u[0];
+
+ assert(!isnan(r[0]));
+ assert(!isnan(r[1]));
+ }
+
+ virtual void write(const size_t step, const Real t, const Real dt, const Tinput& U, const void * const data=nullptr)
+ {
+ m_dumper->write(step, t, dt, U.data(), Tinput::DataType::SIZE*U.size());
+ }
+};
+
+#endif /* RAYLEIGHPLESSET_HBGL_H_EVWDWXSM */
