commit aa2b37100755234f1623e1ee731aeefbe5e82f42
parent 422837ab2f55b7b23bc24297964b5201e17a1815
Author: pvlachas <vlachas@collegium.ethz.ch>
Date: Tue, 9 Feb 2021 14:14:48 +0100
starting with learnedBubbleDynamics in python
Diffstat:
3 files changed, 446 insertions(+), 0 deletions(-)
diff --git a/LearnedBubbleDynamics/config.hjson b/LearnedBubbleDynamics/config.hjson
@@ -0,0 +1,27 @@
+{
+ # time stepper settings
+ ts : "rkv56"
+ ts_tend : 3.0e-6
+ 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 setting
+ N : 4
+ R0 : 10e-6
+ rho:1000
+ nu : 0.0
+ sigma : 0.0
+ gamma : 1.4
+ pv : 0.0
+ p0 : 1.0e5
+ c : 1481.0
+ pamp : 10
+
+ # Initial conditions
+ r0 : 0.01
+ r0_dt : 0.01
+}
diff --git a/LearnedBubbleDynamics/run.py b/LearnedBubbleDynamics/run.py
@@ -0,0 +1,377 @@
+import numpy as np
+import hjson
+import argparse
+
+from timesteper import *
+# Rayleigh-Plesset dynamics (Hilgenfeldt, Brenner, Grossmann, Lohse, 1998)
+
+
+
+
+
+
+class Bubble(object):
+ def __init__(self):
+ self.pos = np.zeros((3))
+
+ def distance(b):
+ dx = b.pos[0] - self.pos[0]
+ dy = b.pos[1] - self.pos[1]
+ dz = b.pos[2] - self.pos[2]
+ return np.sqrt(dx*dx + dy*dy + dz*dz)
+
+class BubbleData(object):
+
+ def __init__(self, params):
+
+ # Number of bubbles and their initial state
+ self.Nbubbles = params["N"]
+
+ self.R0 = np.zeros(self.Nbubbles)
+ self.Rdot0 = np.zeros(self.Nbubbles)
+ self.coords = [Bubble() for i in range(self.Nbubbles)]
+ # np.zeros(self.Nbubbles)
+ self.Dij = np.zeros(self.Nbubbles*self.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;
+
+ for n in range(self.Nbubbles):
+ for i in range(3):
+ self.coords[n].pos[i] = 0.0
+
+ # Gas
+ self.gamma = params["gamma"] # ratio of specific heats for gas
+
+ # Liquid properties
+ self.rhoL = params["rhoL"] # density
+ self.nuL = params["nuL"] # viscosity
+ self.pv = params["pv"] # vaporization pressure
+ self.p0 = params["p0"] # far field pressure
+ self.sigma = params["sigma"] # surface tension
+ self.h = params["h"] # van der Waals hard-core radius
+ self.cl = params["cl"] # speed of sound
+
+ # external forcing
+ self.pamp = params["pamp"]
+ self.frequency = params["frequency"]
+ self.omega = 2.0 * np.pi * self.frequency
+ self.support = params["support"] if params["support"] is not None else 1.0/self.omega
+ self.p_pre = params["p_pre"]
+ self.p_post = params["p_post"]
+ self.t_0 = params["t_0"] if params["t_0"] is not None else 5.0/self.omega
+ self.t_12 = params["t_12"]
+ self.dt_smooth = params["dt_smooth"]
+
+ self.print()
+
+ def pext(self, t):
+ # Driving function
+ pass
+
+ def dpext(self, t):
+ # Driving function
+ pass
+
+ def pBubble(self, U, i):
+ # 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];
+ # pp = (self.p0 - self.pv + 2.0 * sigma / self.R0[i] ) * np.power(self.R0[i] / U[i][0]?, 3.0 * gamma)
+ return pp
+
+ def print(self):
+ # printing bubble data
+ print("Bubble Data:")
+ print("\tNbubbles = {:d}".format(self.Nbubbles))
+ print("\tBubbles:")
+ for i in range(self.Nbubbles):
+ print("\t\tBubble {:d}: R0 = {:e}, Rdot0 = {:e}".format(i, self.R0[i], self.Rdot0[i]))
+ print("\tgamma = {:e}".format(self.gamma))
+ print("\trho = {:e}".format(self.rhoL))
+ print("\tnu = {:e}".format(self.nuL))
+ print("\tpv = {:e}".format(self.pv))
+ print("\tp0 = {:e}".format(self.p0))
+ print("\tsigma = {:e}".format(self.sigma))
+ print("\th = {:e}".format(self.h))
+ print("\tcl = {:e}".format(self.cl))
+ print("\tpamp = {:e}".format(self.pamp))
+ print("\tomega = {:e}".format(self.omega))
+
+
+def getParser():
+ parser = argparse.ArgumentParser()
+ parser.add_argument("--N",
+ help="number of bubbles",
+ type=int,
+ required=False,
+ default=1,
+ )
+ parser.add_argument("--R0",
+ help="initial radius R0 in [m]",
+ type=float,
+ required=False,
+ default=4.0e-6,
+ )
+ parser.add_argument("--Rdot0",
+ help="initial speed in [m/s]",
+ type=float,
+ required=False,
+ default=0.0,
+ )
+
+ # Gas
+ parser.add_argument("--gamma",
+ help="gamma of gas",
+ type=float,
+ required=False,
+ default=1.4,
+ )
+
+ # Liquid
+ parser.add_argument("--rhoL",
+ help="rho of liquid in [kg/m^3]",
+ type=float,
+ required=False,
+ default=1000.0,
+ )
+ parser.add_argument("--nuL",
+ help="nuL of liquid in [m^2/s]",
+ type=float,
+ required=False,
+ default=1.0e-6,
+ )
+ parser.add_argument("--pv",
+ help="pv of liquid in [Pa]",
+ type=float,
+ required=False,
+ default=2400.0,
+ )
+ parser.add_argument("--p0",
+ help="pv of liquid in [Pa]",
+ type=float,
+ required=False,
+ default=1.0e5,
+ )
+ parser.add_argument("--sigma",
+ help="sigma surface tension of liquid in [kg/s^2]",
+ type=float,
+ required=False,
+ default=7.3e-2,
+ )
+ parser.add_argument("--h",
+ help="h of liquid in [kg/s^2]",
+ type=float,
+ required=False,
+ default=4.514673e-7,
+ )
+ parser.add_argument("--cl",
+ help="cl of liquid in [m/s]",
+ type=float,
+ required=False,
+ default=1481.0,
+ )
+
+ # External forcing
+ parser.add_argument("--pamp",
+ help="external pressure amplification",
+ type=float,
+ required=False,
+ default=1.4,
+ )
+ parser.add_argument("--frequency",
+ help="frequency [s^-1]",
+ type=float,
+ required=False,
+ default=26500.0,
+ )
+ parser.add_argument("--support",
+ help="support",
+ type=float,
+ required=False,
+ default=None,
+ )
+ parser.add_argument("--p_pre",
+ help="p_pre",
+ type=float,
+ required=False,
+ default=1.0e5,
+ )
+ parser.add_argument("--p_post",
+ help="p_post",
+ type=float,
+ required=False,
+ default=1.0e5,
+ )
+ parser.add_argument("--t_0",
+ help="t_0",
+ type=float,
+ required=False,
+ default=None,
+ )
+
+ parser.add_argument("--t_12",
+ help="t_12",
+ type=float,
+ required=False,
+ default=1.0e-6,
+ )
+ parser.add_argument("--dt_smooth",
+ help="dt_smooth",
+ type=float,
+ required=False,
+ default=1.0e-9,
+ )
+
+
+
+ parser.add_argument("--ts_aTol",
+ help="ts_aTol",
+ type=float,
+ required=False,
+ default=1.0e-8,
+ )
+ parser.add_argument("--ts_rTol",
+ help="ts_rTol",
+ type=float,
+ required=False,
+ default=1.0e-8,
+ )
+ parser.add_argument("--ts_minScale",
+ help="ts_minScale",
+ type=float,
+ required=False,
+ default=0.2,
+ )
+ parser.add_argument("--ts_maxScale",
+ help="ts_maxScale",
+ type=float,
+ required=False,
+ default=10.0,
+ )
+ parser.add_argument("--ts_alpha",
+ help="ts_alpha",
+ type=float,
+ required=False,
+ default=1.0,
+ )
+ parser.add_argument("--ts_beta",
+ help="ts_beta",
+ type=float,
+ required=False,
+ default=0.0,
+ )
+ parser.add_argument("--ts_safety",
+ help="ts_safety",
+ type=float,
+ required=False,
+ default=0.9,
+ )
+ parser.add_argument("--ts_t0",
+ help="ts_t0",
+ type=float,
+ required=False,
+ default=0.0,
+ )
+ parser.add_argument("--ts_dt",
+ help="ts_dt",
+ type=float,
+ required=False,
+ default=1.0e-4,
+ )
+ parser.add_argument("--ts_nsteps",
+ help="ts_nsteps",
+ type=int,
+ required=False,
+ default=0,
+ )
+ parser.add_argument("--ts_tend",
+ help="ts_tend",
+ type=float,
+ required=False,
+ default=0,
+ )
+ parser.add_argument("--ts_dtDump",
+ help="ts_dtDump",
+ type=float,
+ required=False,
+ default=-1.0,
+ )
+ parser.add_argument("--ts_writeGranularity",
+ help="ts_writeGranularity",
+ type=int,
+ required=False,
+ default=1,
+ )
+ parser.add_argument("--ts_reportGranularity",
+ help="ts_reportGranularity",
+ type=int,
+ required=False,
+ default=100,
+ )
+ parser.add_argument("--ts_restart",
+ help="ts_restart",
+ type=bool,
+ required=False,
+ default=False,
+ )
+ parser.add_argument("--ts_restart_step",
+ help="ts_restart_step",
+ type=int,
+ required=False,
+ default=0,
+ )
+
+
+ return parser
+
+
+
+class RayleighPlesset(object):
+
+ def compute(U, rhs, t, bd):
+
+ r[0] = u[1]
+
+ u0_power_3 = np.power(u[0], 3)
+ dR03 = np.power(bd.R0, 3) - np.power(bd.h, 3)
+ dR3 = u0_power_3 - np.power(bd.h, 3)
+
+ r[1] = bd.rInv * (bd.pv - bd.p0 - bd.pext(t))
+ r[1] += bd.rInv*(bd.p0 - bd.pv + 2.0*bd.sigma/bd.R0[0])* np.power(dR03/dR3,bd.gamma)*(1.0 - 3.0*bd.gamma*u0_power_3*u[1]/(bd.cl*dR3))
+ r[1] -= 1.5 * (u[1] * u[1])
+ r[1] -= 4.0 * bd.nuL * u[1]/u[0]
+ r[1] -= 2.0 * bd.sigma * bd.rInv / u[0]
+ r[1] /= u[0]
+
+
+
+
+def main():
+ parser = getParser()
+ args = parser.parse_args()
+ args_default_dict = args.__dict__
+ args_json_dict = hjson.load(open("config.hjson", "rb"))
+ args_dict = args_default_dict.copy()
+ args_dict.update(args_json_dict)
+
+ bd = BubbleData(args_dict)
+
+ time_stepper = Stepper(args_dict)
+
+
+
+if __name__ == '__main__':
+ main()
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/LearnedBubbleDynamics/timesteper.py b/LearnedBubbleDynamics/timesteper.py
@@ -0,0 +1,42 @@
+
+
+class Stepper(object):
+ def __init__(self, params):
+ # Stepper properties
+ self.aTol = params["ts_aTol"] # absolute tolerance
+ self.rTol = params["ts_rTol"] # relative tolerance
+
+ self.minScale = params["ts_minScale"]
+ self.maxScale = params["ts_maxScale"]
+ self.alpha = params["ts_alpha"]
+ self.beta = params["ts_beta"]
+ self.safety = params["ts_safety"]
+ self.t = params["ts_t0"]
+ self.dt = params["ts_dt"]
+ self.step = 0
+ self.nsteps = params["ts_nsteps"]
+ self.tFinal = params["ts_tend"]
+ self.dtDump = params["ts_dtDump"]
+ self.tDump = params["pv"]
+ self.writeGranularity = params["ts_writeGranularity"]
+ self.writeCount = 0
+ self.reportGranularity = params["ts_reportGranularity"]
+ self.bFixedStep = True
+ # Restarts
+ self.bRestart = params["ts_restart"]
+ self.restartstep = params["ts_restart_step"]
+ self.print()
+
+ def print(self):
+ print("Time Stepper :")
+ print("\tAbsolute tolerance = {:e}".format(self.aTol))
+ print("\tRelative tolerance = {:e}".format(self.rTol))
+ print("\tMinimum time step scale = {:e}".format(self.minScale))
+ print("\tMaximum time step scale = {:e}".format(self.maxScale))
+ print("\tError PI control alpha = {:e}".format(self.alpha))
+ print("\tError PI control beta = {:e}".format(self.beta))
+
+def getParser():
+ parser = argparse.ArgumentParser()
+
+
