From 672e26c0f6625062e37bed2743dafff0277b488a Mon Sep 17 00:00:00 2001 From: Alice Harpole Date: Tue, 9 Jul 2019 15:24:29 -0400 Subject: [PATCH] sph solver --- pyro.py | 1 + sph/__init__.py | 9 + sph/_defaults | 12 + sph/compute.py | 52 ++++ sph/particles.py | 570 ++++++++++++++++++++++++++++++++++ sph/problems/__init__.py | 1 + sph/problems/_circle.defaults | 1 + sph/problems/circle.py | 59 ++++ sph/problems/inputs.circle | 31 ++ sph/simulation.py | 313 +++++++++++++++++++ 10 files changed, 1049 insertions(+) create mode 100644 sph/__init__.py create mode 100644 sph/_defaults create mode 100644 sph/compute.py create mode 100644 sph/particles.py create mode 100644 sph/problems/__init__.py create mode 100644 sph/problems/_circle.defaults create mode 100644 sph/problems/circle.py create mode 100644 sph/problems/inputs.circle create mode 100644 sph/simulation.py diff --git a/pyro.py b/pyro.py index 67de7b688..9d5d50d4e 100755 --- a/pyro.py +++ b/pyro.py @@ -24,6 +24,7 @@ "diffusion", "incompressible", "lm_atm", + "sph", "swe"] diff --git a/sph/__init__.py b/sph/__init__.py new file mode 100644 index 000000000..986b63a98 --- /dev/null +++ b/sph/__init__.py @@ -0,0 +1,9 @@ +""" +The pyro swe hydrodynamics solver. This implements the +second-order (piecewise-linear), unsplit method of Colella 1990. + +""" + +__all__ = ["simulation"] + +from .simulation import Simulation diff --git a/sph/_defaults b/sph/_defaults new file mode 100644 index 000000000..201aacd6a --- /dev/null +++ b/sph/_defaults @@ -0,0 +1,12 @@ +[driver] +cfl = 0.8 + + +[sph] +np = 100 +grav = -9.8 +particle_size = 5.e-2 +rho0 = 1.e3 +k = 1.e3 +mu = 0.1 +dt = 1.e-4 diff --git a/sph/compute.py b/sph/compute.py new file mode 100644 index 000000000..6d1fa56c3 --- /dev/null +++ b/sph/compute.py @@ -0,0 +1,52 @@ +import numpy as np + + +def compute_density(my_data, ivars): + U = my_data.data + h = my_data.get_aux("h") + + U[:, ivars.irho] = 0 + + for i in range(my_data.np): + + C = 4 / (np.pi * h**8) + U[i, ivars.irho] += 4 * U[i, ivars.im] / (np.pi * h**2) + + for j in range(i+1): + r_ij = U[i, ivars.ix:ivars.iy + 1] - U[j, ivars.ix:ivars.iy + 1] + z = h**2 - np.dot(r_ij, r_ij) + + if z > 0: + U[i, ivars.irho] += C * U[i, ivars.im] * z**3 + U[j, ivars.irho] += C * U[j, ivars.im] * z**3 + +def compute_acceleration(my_data, ivars): + + compute_density(my_data, ivars) + + U = my_data.data + h = my_data.get_aux("h") + rho0 = my_data.get_aux("rho0") + mu = my_data.get_aux("mu") + k = my_data.get_aux("k") + + U[:, ivars.iax] = 0 + U[:, ivars.iay] = my_data.get_aux("g") + + for i in range(my_data.np): + for j in range(i): + r_ij = U[i, ivars.ix:ivars.iy + 1] - U[j, ivars.ix:ivars.iy + 1] + v_ij = U[i, ivars.iu:ivars.iv + 1] - U[j, ivars.iu:ivars.iv + 1] + q = np.sqrt(np.dot(r_ij, r_ij)) / h + if q > 0: + f_ij = (1 - q) * (15 * k * (U[i, ivars.irho] + + U[j, ivars.irho] - 2 * rho0) * + (1 - q) / q * r_ij - 40 * mu * v_ij) + else: + f_ij = - 40 * mu * v_ij + + U[i, ivars.iax:ivars.iay + 1] += f_ij * U[j, ivars.im] / \ + (np.pi * h**4 * U[j, ivars.irho] * U[i, ivars.irho]) + + U[j, ivars.iax:ivars.iay + 1] -= f_ij * U[i, ivars.im] / \ + (np.pi * h**4 * U[j, ivars.irho] * U[i, ivars.irho]) diff --git a/sph/particles.py b/sph/particles.py new file mode 100644 index 000000000..1e424c018 --- /dev/null +++ b/sph/particles.py @@ -0,0 +1,570 @@ +""" +The patch module defines the classes necessary to describe finite-volume +data and the grid that it lives on. + +Typical usage: + +* create the grid:: + + grid = Domain2d(nx, ny) + +* create the data that lives on that grid:: + + data = ParticleData2d(grid) + + bc = BC(xlb="reflect", xrb="reflect", + ylb="outflow", yrb="outflow") + data.register_var("density", bc) + ... + + data.create() + +* initialize some data:: + + dens = data.get_var("density") + dens[:, :] = ... + + +* fill the ghost cells:: + + data.fill_BC("density") + +""" +from __future__ import print_function + +import numpy as np + +import h5py + +from util import msg + +from sph.compute import compute_density + +class Variables(object): + """ + a container class for easy access to the different sph + variables by an integer key + """ + def __init__(self, myd): + self.nvar = len(myd.names) + + # conserved variables -- we set these when we initialize for + # they match the CellCenterData2d object + self.im = myd.names.index("mass") + self.irho = myd.names.index("density") + self.ix = myd.names.index("x-position") + self.iy = myd.names.index("y-position") + self.iu = myd.names.index("x-velocity") + self.iv = myd.names.index("y-velocity") + self.iuh = myd.names.index("half-x-velocity") + self.ivh = myd.names.index("half-y-velocity") + self.iax = myd.names.index("x-acceleration") + self.iay = myd.names.index("y-acceleration") + +class Domain2d(object): + """ + the 2-d grid class. The grid object will contain the coordinate + information (at various centerings). + + A basic (1-d) representation of the layout is:: + + | | | X | | | | X | | | + +--*--+- // -+--*--X--*--+--*--+- // -+--*--+--*--X--*--+- // -+--*--+ + 0 ng-1 ng ng+1 ... ng+nx-1 ng+nx 2ng+nx-1 + + ilo ihi + + |<- ng guardcells->|<---- nx interior zones ----->|<- ng guardcells->| + + The '*' marks the data locations. + """ + + # pylint: disable=too-many-instance-attributes + + def __init__(self, xmin=0.0, xmax=1.0, ymin=0.0, ymax=1.0): + """ + Create a Domain2d object. + + The only data that we require is the number of points that + make up the mesh in each direction. Optionally we take the + extrema of the domain (default is [0,1]x[0,1]) and number of + ghost cells (default is 1). + + Note that the Domain2d object only defines the discretization, + it does not know about the boundary conditions, as these can + vary depending on the variable. + + Parameters + ---------- + xmin : float, optional + Physical coordinate at the lower x boundary + xmax : float, optional + Physical coordinate at the upper x boundary + ymin : float, optional + Physical coordinate at the lower y boundary + ymax : float, optional + Physical coordinate at the upper y boundary + """ + + # domain extrema + self.xmin = xmin + self.xmax = xmax + + self.ymin = ymin + self.ymax = ymax + + def __eq__(self, other): + """ are two domains equivalent? """ + result = (self.xmin == other.xmin and self.xmax == other.xmax and + self.ymin == other.ymin and self.ymax == other.ymax) + + return result + +class ParticleData2d(object): + """ + A class to define cell-centered data that lives on a grid. A + ParticleData2d object is built in a multi-step process before + it can be used. + + * Create the object. We pass in a grid object to describe where + the data lives:: + + my_data = patch.ParticleData2d(myGrid) + + * Register any variables that we expect to live on this patch. + Here BC describes the boundary conditions for that variable:: + + my_data.register_var('density', BC) + my_data.register_var('x-momentum', BC) + ... + + * Register any auxillary data -- these are any parameters that are + needed to interpret the data outside of the simulation (for + example, the gamma for the equation of state):: + + my_data.set_aux(keyword, value) + + * Finish the initialization of the patch:: + + my_data.create() + + This last step actually allocates the storage for the state + variables. Once this is done, the patch is considered to be + locked. New variables cannot be added. + """ + + # pylint: disable=too-many-instance-attributes + + def __init__(self, domain, rp, bc, dtype=np.float64): + + """ + Initialize the ParticleData2d object. + + Parameters + ---------- + domain : Domain2d object + The domain upon which the data will live + dtype : NumPy data type, optional + The datatype of the data we wish to create (defaults to + np.float64 + """ + + self.domain = domain + self.np = rp.get_param("sph.np") + + self.dtype = dtype + self.data = None + + self.names = [] + self.vars = self.names # backwards compatibility hack + self.nvar = 0 + + self.aux = {} + + self.BCs = bc + + # time + self.t = -1.0 + + self.initialized = 0 + + def register_var(self, name): + """ + Register a variable with ParticleData2d object. + + Parameters + ---------- + name : str + The variable name + bc : BC object + The boundary conditions that describe the actions to take + for this variable at the physical domain boundaries. + """ + + if self.initialized == 1: + msg.fail("ERROR: domain already initialized") + + self.names.append(name) + self.nvar += 1 + + def set_aux(self, keyword, value): + """ + Set any auxillary (scalar) data. This data is simply carried + along with the ParticleData2d object + + Parameters + ---------- + keyword : str + The name of the datum + value : any time + The value to associate with the keyword + """ + self.aux[keyword] = value + + def create(self): + """ + Called after all the variables are registered and allocates + the storage for the state data. + """ + + if self.initialized == 1: + msg.fail("ERROR: domain already initialized") + + self.data = np.zeros((self.np, self.nvar)) + + self.initialized = 1 + + def __str__(self): + """ print out some basic information about the ParticleData2d + object """ + + if self.initialized == 0: + my_str = "ParticleData2d object not yet initialized" + return my_str + + my_str = " nvars = {}\n".format(self.nvar) + my_str += " variables:\n" + + for n in range(self.nvar): + my_str += "%16s: min: %15.10f max: %15.10f\n" % \ + (self.names[n], self.min(self.names[n]), self.max(self.names[n])) + # my_str += "%16s BCs: -x: %-12s +x: %-12s -y: %-12s +y: %-12s\n" %\ + # (" ", self.BCs[self.names[n]].xlb, + # self.BCs[self.names[n]].xrb, + # self.BCs[self.names[n]].ylb, + # self.BCs[self.names[n]].yrb) + + return my_str + + def get_var(self, name): + """ + Return a data array for the variable described by name. Stored + variables will be checked first, and then any derived variables + will be checked. + + For a stored variable, changes made to this are automatically + reflected in the ParticleData2d object. + + Parameters + ---------- + name : str + The name of the variable to access + + Returns + ------- + out : ndarray + The array of data corresponding to the variable name + + """ + try: + n = self.names.index(name) + except ValueError: + for f in self.derives: + var = f(self, name) + if len(var) > 0: + return var + raise KeyError("name {} is not valid".format(name)) + else: + return self.data[:, n] + + def get_var_by_index(self, n): + """ + Return a data array for the variable with index n in the + data array. Any changes made to this are automatically + reflected in the ParticleData2d object. + + Parameters + ---------- + n : int + The index of the variable to access + + Returns + ------- + out : ndarray + The array of data corresponding to the index + + """ + return self.data[:, n] + + def get_vars(self): + """ + Return the entire data array. Any changes made to this + are automatically reflected in the ParticleData2d object. + + Returns + ------- + out : ndarray + The array of data + + """ + return self.data + + def get_aux(self, keyword): + """ + Get the auxillary data associated with keyword + + Parameters + ---------- + keyword : str + The name of the auxillary data to access + + Returns + ------- + out : variable type + The value corresponding to the keyword + + """ + if keyword in self.aux.keys(): + return self.aux[keyword] + + return None + + def zero(self, name): + """ + Zero out the data array associated with variable name. + + Parameters + ---------- + name : str + The name of the variable to zero + + """ + n = self.names.index(name) + self.data[:, n] = 0.0 + + def fill_BC_all(self): + """ + Fill boundary conditions on all variables. + """ + ivars = Variables(self) + + bc = self.BCs + U = self.data + + # -x boundary + if bc.xlb in ["outflow", "neumann"]: + # in the future we could just kill these particles + pass + elif bc.xlb in ["reflect-even", "reflect-odd"]: + # print("hi i'm lower x reflective") + # tbounce = (U[:, ivars.ix] - self.domain.xmin) / U[:, ivars.iu] + for p in range(self.np): + if U[p, ivars.ix] < self.domain.xmin: + U[p, ivars.ix] = 2 * self.domain.xmin - U[p, ivars.ix] + U[p, ivars.iu] = -U[p, ivars.iu] + U[p, ivars.iuh] = -U[p, ivars.iuh] + elif bc.xlb == "periodic": + # print("hi i'm lower x periodic") + for p in range(self.np): + if U[p, ivars.ix] < self.domain.xmin: + U[p, ivars.ix] = self.domain.xmax - (self.domain.xmin - U[p, ivars.ix]) + + # +x boundary + if bc.xrb in ["outflow", "neumann"]: + # in the future we could just kill these particles + pass + elif bc.xrb in ["reflect-even", "reflect-odd"]: + for p in range(self.np): + if U[p, ivars.ix] > self.domain.xmax: + U[p, ivars.ix] = 2 * self.domain.xmax - U[p, ivars.ix] + U[p, ivars.iu] = -U[p, ivars.iu] + U[p, ivars.iuh] = -U[p, ivars.iuh] + elif bc.xrb == "periodic": + for p in range(self.np): + if U[p, ivars.ix] > self.domain.xmax: + U[p, ivars.ix] = self.domain.xmin + (U[p, ivars.ix] - self.domain.xmax) + + # -y boundary + if bc.ylb in ["outflow", "neumann"]: + # in the future we could just kill these particles + pass + elif bc.ylb in ["reflect-even", "reflect-odd"]: + for p in range(self.np): + if U[p, ivars.iy] < self.domain.ymin: + U[p, ivars.iy] = 2 * self.domain.ymin - U[p, ivars.iy] + U[p, ivars.iv] = -U[p, ivars.iv] + U[p, ivars.ivh] = -U[p, ivars.ivh] + elif bc.ylb == "periodic": + for p in range(self.np): + if U[p, ivars.iy] < self.domain.ymin: + U[p, ivars.iy] = self.domain.ymax - (self.domain.ymin - U[p, ivars.iy]) + + # +y boundary + if bc.yrb in ["outflow", "neumann"]: + # in the future we could just kill these particles + pass + elif bc.yrb in ["reflect-even", "reflect-odd"]: + for p in range(self.np): + if U[p, ivars.iy] > self.domain.ymax: + U[p, ivars.iy] = 2 * self.domain.ymax - U[p, ivars.iy] + U[p, ivars.iv] = -U[p, ivars.iv] + U[p, ivars.ivh] = -U[p, ivars.ivh] + elif bc.yrb == "periodic": + for p in range(self.np): + if U[p, ivars.iy] > self.domain.ymax: + U[p, ivars.iy] = self.domain.ymin + (U[p, ivars.iy] - self.domain.ymax) + + + def min(self, name): + """ + return the minimum of the variable name in the domain's valid region + """ + n = self.names.index(name) + return np.min(self.data[:, n]) + + def max(self, name): + """ + return the maximum of the variable name in the domain's valid region + """ + n = self.names.index(name) + return np.max(self.data[:, n]) + + def write(self, filename): + """ + create an output file in HDF5 format and write out our data and + grid. + """ + + if not filename.endswith(".h5"): + filename += ".h5" + + with h5py.File(filename, "w") as f: + self.write_data(f) + + def write_data(self, f): + """ + write the data out to an hdf5 file -- here, f is an h5py + File pbject + + """ + + # auxillary data + gaux = f.create_group("aux") + for k, v in self.aux.items(): + gaux.attrs[k] = v + + # grid information + ggrid = f.create_group("grid") + + # data + gstate = f.create_group("state") + + # for n in range(self.nvar): + # gvar = gstate.create_group(self.names[n]) + # gvar.create_dataset("data", + # data=self.get_var_by_index(n).v()) + # gvar.attrs["xlb"] = self.BCs[self.names[n]].xlb + # gvar.attrs["xrb"] = self.BCs[self.names[n]].xrb + # gvar.attrs["ylb"] = self.BCs[self.names[n]].ylb + # gvar.attrs["yrb"] = self.BCs[self.names[n]].yrb + + def pretty_print(self, var, fmt=None): + """print out the contents of the data array with pretty formatting + indicating where ghost cells are.""" + a = self.get_var(var) + a.pretty_print(fmt=fmt) + + def normalize_mass(self): + ivars = Variables(self) + + U = self.data + + compute_density(self, ivars) + + rho0 = self.get_aux("rho0") + rho2s = np.sum(U[:, ivars.irho]**2) + rhos = np.sum(U[:, ivars.irho]) + + U[:, ivars.im] *= rho0 * rhos / rho2s + +def particle_data_clone(old): + """ + Create a new ParticleData2d object that is a copy of an existing + one + + Parameters + ---------- + old : ParticleData2d object + The ParticleData2d object we wish to copy + + Note + ---- + It may be that this whole thing can be replaced with a copy.deepcopy() + + """ + + if not isinstance(old, ParticleData2d): + msg.fail("Can't clone object") + + # we may be a type derived from ParticleData2d, so use the same + # type + myt = type(old) + new = myt(old.domain, old.np, dtype=old.dtype) + + for n in range(old.nvar): + new.register_var(old.names[n], old.BCs[old.names[n]]) + + new.create() + + new.aux = old.aux.copy() + new.data = old.data.copy() + + return new + + +def do_demo(): + """ show examples of the patch methods / classes """ + + pass + + # import util.io as io + # + # # illustrate basic mesh operations + # + # myg = Domain2d(8, 16, xmax=1.0, ymax=2.0) + # + # mydata = ParticleData2d(myg) + # + # bc = bnd.BC() + # + # mydata.register_var("a", bc) + # mydata.create() + # + # a = mydata.get_var("a") + # a[:, :] = np.exp(-(myg.x2d - 0.5)**2 - (myg.y2d - 1.0)**2) + # + # print(mydata) + # + # # output + # print("writing\n") + # mydata.write("mesh_test") + # + # print("reading\n") + # myd2 = io.read("mesh_test") + # print(myd2) + # + # mydata.pretty_print("a") + + +if __name__ == "__main__": + do_demo() diff --git a/sph/problems/__init__.py b/sph/problems/__init__.py new file mode 100644 index 000000000..c0c82ebdd --- /dev/null +++ b/sph/problems/__init__.py @@ -0,0 +1 @@ +__all__ = ['acoustic_pulse', 'advect', 'dam', 'kh', 'logo', 'quad', 'test'] diff --git a/sph/problems/_circle.defaults b/sph/problems/_circle.defaults new file mode 100644 index 000000000..57ac7ee0a --- /dev/null +++ b/sph/problems/_circle.defaults @@ -0,0 +1 @@ +[circle] diff --git a/sph/problems/circle.py b/sph/problems/circle.py new file mode 100644 index 000000000..4c8090962 --- /dev/null +++ b/sph/problems/circle.py @@ -0,0 +1,59 @@ +from __future__ import print_function + +import sys +import numpy as np +import mesh.patch as patch +from util import msg +from sph.particles import Variables, ParticleData2d + + +def init_data(myd, rp): + """initialize the circle problem.""" + + msg.bold("initializing the circle problem...") + + # make sure that we are passed a valid patch object + if not isinstance(myd, ParticleData2d): + print("ERROR: patch invalid in circle.py") + print(myd.__class__) + sys.exit() + + U = myd.data + ivars = Variables(myd) + + h = myd.get_aux("h") + hh = h / 1.8 + + count = 0 + + def circ_indicator(x, y): + dx = x - 0.5 + dy = y - 0.5 + + return dx**2 + dy**2 < 0.25**2 + + for x in np.arange(0, 1, hh): + for y in np.arange(0, 1, hh): + count += int(circ_indicator(x, y)) + + p = 0 + + for x in np.arange(0, 1, hh): + for y in np.arange(0, 1, hh): + if circ_indicator(x, y): + U[p, ivars.ix] = x + U[p, ivars.iy] = y + U[p, ivars.iu:ivars.iv + 1] = 0 + p += 1 + + if len(U[:, 0] > p): + myd.np = p + myd.data = U[:p, :] + + # give all particles the same mass + U[:, ivars.im] = 1 + + +def finalize(): + """ print out any information to the user at the end of the run """ + pass diff --git a/sph/problems/inputs.circle b/sph/problems/inputs.circle new file mode 100644 index 000000000..c44eecdcc --- /dev/null +++ b/sph/problems/inputs.circle @@ -0,0 +1,31 @@ +[driver] +max_steps = 200 +tmax = 0.3 + +[sph] +np = 500 +grav = -9.8 +particle_size = 5.e-2 +rho0 = 1.e3 +k = 1.e3 +mu = 0.1 +dt = 1.e-4 + + +[io] +basename = circle_ + + +[mesh] +xmax = 1.0 +ymax = 1.0 + +xlboundary = periodic +xrboundary = periodic + +ylboundary = reflect +yrboundary = reflect + + +[vis] +dovis = 1 diff --git a/sph/simulation.py b/sph/simulation.py new file mode 100644 index 000000000..b9ab0008f --- /dev/null +++ b/sph/simulation.py @@ -0,0 +1,313 @@ +from __future__ import print_function + +import importlib + +import numpy as np +import matplotlib.pyplot as plt + +import sph.compute as compute +import mesh.boundary as bnd +from simulation_null import NullSimulation +import util.plot_tools as plot_tools +from util import msg +import sph.particles as particles +from scipy.interpolate import interp2d + +def domain_setup(rp, ng=1): + + try: + xmin = rp.get_param("mesh.xmin") + except KeyError: + xmin = 0.0 + msg.warning("mesh.xmin not set, defaulting to 0.0") + + try: + xmax = rp.get_param("mesh.xmax") + except KeyError: + xmax = 1.0 + msg.warning("mesh.xmax not set, defaulting to 1.0") + + try: + ymin = rp.get_param("mesh.ymin") + except KeyError: + ymin = 0.0 + msg.warning("mesh.ymin not set, defaulting to 0.0") + + try: + ymax = rp.get_param("mesh.ymax") + except KeyError: + ymax = 1.0 + msg.warning("mesh.ynax not set, defaulting to 1.0") + + my_domain = particles.Domain2d(xmin=xmin, xmax=xmax, + ymin=ymin, ymax=ymax) + return my_domain + +def bc_setup(rp): + + # first figure out the BCs + try: + xlb_type = rp.get_param("mesh.xlboundary") + except KeyError: + xlb_type = "periodic" + msg.warning("mesh.xlboundary is not set, defaulting to periodic") + + try: + xrb_type = rp.get_param("mesh.xrboundary") + except KeyError: + xrb_type = "periodic" + msg.warning("mesh.xrboundary is not set, defaulting to periodic") + + try: + ylb_type = rp.get_param("mesh.ylboundary") + except KeyError: + ylb_type = "periodic" + msg.warning("mesh.ylboundary is not set, defaulting to periodic") + + try: + yrb_type = rp.get_param("mesh.yrboundary") + except KeyError: + yrb_type = "periodic" + msg.warning("mesh.yrboundary is not set, defaulting to periodic") + + bc = bnd.BC(xlb=xlb_type, xrb=xrb_type, + ylb=ylb_type, yrb=yrb_type) + + return bc + + +class Simulation(NullSimulation): + """The main simulation class for the corner transport upwind + sph hydrodynamics solver + + """ + + def __init__(self, solver_name, problem_name, rp, timers=None, data_class=particles.ParticleData2d): + + return super().__init__(solver_name, problem_name, rp, timers, data_class) + + def initialize(self, extra_vars=None, ng=4): + """ + Initialize the grid and variables for sph flow and set + the initial conditions for the chosen problem. + """ + my_domain = domain_setup(self.rp) + bc = bc_setup(self.rp) + + my_data = self.data_class(my_domain, self.rp, bc) + + # are we dealing with solid boundaries? we'll use these for + # the Riemann solver + self.solid = bnd.bc_is_solid(bc) + + # density and energy + my_data.register_var("mass") + my_data.register_var("density") + my_data.register_var("x-position") + my_data.register_var("y-position") + my_data.register_var("x-velocity") + my_data.register_var("y-velocity") + my_data.register_var("half-x-velocity") + my_data.register_var("half-y-velocity") + my_data.register_var("x-acceleration") + my_data.register_var("y-acceleration") + + # store the gravitational acceration g as an auxillary quantity + # so we can have a + # self-contained object stored in output files to make plots. + # store grav because we'll need that in some BCs + my_data.set_aux("g", self.rp.get_param("sph.grav")) + my_data.set_aux("h", self.rp.get_param("sph.particle_size")) + my_data.set_aux("rho0", self.rp.get_param("sph.rho0")) + my_data.set_aux("k", self.rp.get_param("sph.k")) + my_data.set_aux("mu", self.rp.get_param("sph.mu")) + + my_data.create() + + self.cc_data = my_data + + self.ivars = particles.Variables(my_data) + + # initial conditions for the problem + problem = importlib.import_module("{}.problems.{}".format( + self.solver_name, self.problem_name)) + problem.init_data(self.cc_data, self.rp) + + self.cc_data.normalize_mass() + + self.compute_timestep() + + if self.verbose > 0: + print(my_data) + + def compute_timestep(self): + """ + The timestep function computes the advective timestep (CFL) + constraint. The CFL constraint says that information cannot + propagate further than one zone per timestep. + + We use the driver.cfl parameter to control what fraction of the + CFL step we actually take. + """ + + cfl = self.rp.get_param("driver.cfl") + + dx = self.cc_data.domain.xmax - self.cc_data.domain.xmin + dy = self.cc_data.domain.ymax - self.cc_data.domain.ymin + + max_u = np.max(np.abs(self.cc_data.get_var("x-velocity"))) + max_v = np.max(np.abs(self.cc_data.get_var("y-velocity"))) + + self.dt = np.min((self.rp.get_param("sph.dt"), cfl*dx/max_u, cfl*dy/max_v)) + + def preevolve(self): + """ + Evolve the equations of sph hydrodynamics through a + timestep dt. + """ + + tm_evolve = self.tc.timer("preevolve") + tm_evolve.begin() + + U = self.cc_data.data + ivars = particles.Variables(self.cc_data) + + # compute the acceleration + compute.compute_acceleration(self.cc_data, ivars) + + # update velocity at half timestep + U[:, ivars.iuh:ivars.ivh+1] = U[:, ivars.iu:ivars.iv+1] + U[:, ivars.iax:ivars.iay+1] * self.dt * 0.5 + + # update velocity + U[:, ivars.iu:ivars.iv+1] += U[:, ivars.iax:ivars.iay+1] * self.dt + + # update position + U[:, ivars.ix:ivars.iy+1] += U[:, ivars.iuh:ivars.ivh+1] * self.dt + + tm_evolve.end() + + self.cc_data.fill_BC_all() + + def evolve(self): + """ + Evolve the equations of sph hydrodynamics through a + timestep dt. + """ + + tm_evolve = self.tc.timer("evolve") + tm_evolve.begin() + + U = self.cc_data.data + ivars = particles.Variables(self.cc_data) + + # compute the acceleration + compute.compute_acceleration(self.cc_data, ivars) + + # update velocity at half timestep + U[:, ivars.iuh:ivars.ivh+1] += U[:, ivars.iax:ivars.iay+1] * self.dt + + # update velocity + U[:, ivars.iu:ivars.iv+1] = U[:, ivars.iuh:ivars.ivh+1] + U[:, ivars.iax:ivars.iay+1] * self.dt * 0.5 + + # update position + U[:, ivars.ix:ivars.iy+1] += U[:, ivars.iuh:ivars.ivh+1] * self.dt + + # increment the time + self.cc_data.t += self.dt + self.n += 1 + + self.cc_data.fill_BC_all() + + tm_evolve.end() + + def dovis(self): + """ + Do runtime visualization. + """ + + plt.clf() + + plt.rc("font", size=10) + + _, axes, cbar_title = plot_tools.setup_axes(self.cc_data.domain, 1) + + myd = self.cc_data + # ivars = particles.Variables(self.cc_data) + + xs = myd.get_var("x-position") + ys = myd.get_var("y-position") + + # compute.compute_density(self.cc_data, ivars) + # + # dens = myd.get_var("density") + # + # # f = interp2d(xs, ys, dens, kind='cubic') + # + # X = np.linspace(self.cc_data.domain.xmin, self.cc_data.domain.xmax, 100) + # Y = np.linspace(self.cc_data.domain.ymin, self.cc_data.domain.ymax, 100) + + # interp_dens = f(X, Y) + + # print(interp_dens) + + ax = axes[0] + + # ax.imshow(np.transpose(interp_dens)) + ax.scatter(xs, ys)#, c=np.array(range(self.cc_data.np))) + ax.set_xlim([self.cc_data.domain.xmin, self.cc_data.domain.xmax]) + ax.set_ylim([self.cc_data.domain.ymin, self.cc_data.domain.ymax]) + + ax.set_xlabel("x") + ax.set_ylabel("y") + + # # access g from the cc_data object so we can use dovis + # # outside of a running simulation. + # g = self.cc_data.get_aux("g") + # + # + # h = q[:, :, ivars.ih] + # u = q[:, :, ivars.iu] + # v = q[:, :, ivars.iv] + # fuel = q[:, :, ivars.ix] + # + # magvel = np.sqrt(u**2 + v**2) + # + # myg = self.cc_data.grid + # + # vort = myg.scratch_array() + # + # dv = 0.5*(v.ip(1) - v.ip(-1))/myg.dx + # du = 0.5*(u.jp(1) - u.jp(-1))/myg.dy + # + # vort.v()[:, :] = dv - du + # + # fields = [h, magvel, fuel, vort] + # field_names = [r"$h$", r"$|U|$", r"$X$", r"$\nabla\times U$"] + # + # _, axes, cbar_title = plot_tools.setup_axes(myg, len(fields)) + # + # for n, ax in enumerate(axes): + # v = fields[n] + # + # img = ax.imshow(np.transpose(v.v()), + # interpolation="nearest", origin="lower", + # extent=[myg.xmin, myg.xmax, myg.ymin, myg.ymax], + # cmap=self.cm) + # + # ax.set_xlabel("x") + # ax.set_ylabel("y") + # + # # needed for PDF rendering + # cb = axes.cbar_axes[n].colorbar(img) + # cb.solids.set_rasterized(True) + # cb.solids.set_edgecolor("face") + # + # if cbar_title: + # cb.ax.set_title(field_names[n]) + # else: + # ax.set_title(field_names[n]) + # + plt.figtext(0.05, 0.0125, "t = {:10.5g}".format(self.cc_data.t)) + + plt.pause(0.0001) + plt.draw()