A tour of a Spheral script

Running a problem with Spheral involves assembling the necessary objects in Python: typically a set of physics packages, equations of state, materials, time integrator, etc. These can then be handed to an object designed to coordinate and run the problem (the SpheralController), generating output such as visualization files or invoking custom analysis routines the user might provide. Each of these concepts will be discussed in depth later in this manual, but first let’s take a tour of a simple example designed to run a classic hydrodynamics test: the Taylor-Sedov blastwave [1].

The analytical description of this problem starts with an infinitesimal energy spike deposited in a uniform pressureless gas. This results in a shock or blastwave that propagates radially away from the original energy spike, and can be described by a self-similar analytic solution for a simple analytical material (such as a \(\gamma\) law gas).

../../../_images/Sedov_demo.gif

Time sequence of the mass density evolution in the Sedov problem using default Spheral visualization via Visit. The black dots show the positions of the SPH nodes, while the spatial color map of density is produced by coloring polygonal Voronoi cells containing and constructed from the node positions.

First let’s look at the example Python script in its entirety, and then we’ll go through it section by section in more detail:

  1#-------------------------------------------------------------------------------
  2# The cylindrical Sedov test case (2-D).
  3#-------------------------------------------------------------------------------
  4import os, shutil
  5import mpi
  6from Spheral2d import *
  7from SpheralTestUtilities import *
  8from GenerateNodeDistribution2d import *
  9from PeanoHilbertDistributeNodes import distributeNodes2d
 10from SpheralMatplotlib import *
 11
 12title("2-D SPH hydrodynamics demonstration of the cylindrical Sedov problem")
 13
 14#-------------------------------------------------------------------------------
 15# Generic problem parameters
 16#-------------------------------------------------------------------------------
 17commandLine(
 18    # Geometry of the problem
 19    rmin = 0.0,
 20    rmax = 1.0,
 21
 22    # Resolution and how to lay down the points
 23    nRadial = 50,
 24    nPerh = 4.01,
 25
 26    # Initial conditions
 27    rho0 = 1.0,
 28    eps0 = 0.0,
 29    Espike = 0.25*1.0,
 30
 31    # Material equation of state options (gamma-law gas in this case)
 32    gamma = 5.0/3.0,
 33    mu = 1.0,
 34
 35    # Time stepping
 36    goalTime = 0.5,
 37    dt = 1.0e-8,
 38    dtGrowth = 2.0,
 39
 40    # IO
 41    vizCycle = None,
 42    vizTime = 0.1,
 43    restoreCycle = -1,
 44    restartStep = 1000,
 45    clearDirectories = False,
 46    dataDirBase = "dumps-cylindrical-Sedov",
 47    graphics = True,
 48)
 49
 50#-------------------------------------------------------------------------------
 51# Path names.
 52#-------------------------------------------------------------------------------
 53dataDir = os.path.join(dataDirBase,
 54                       "nr={}".format(nRadial))
 55restartDir = os.path.join(dataDir, "restarts")
 56vizDir = os.path.join(dataDir, "visit")
 57restartBaseName = os.path.join(restartDir, "Sedov-cylindrical-2d-%i" % nRadial)
 58
 59#-------------------------------------------------------------------------------
 60# Check if the necessary output directories exist.  If not, create them.
 61#-------------------------------------------------------------------------------
 62if mpi.rank == 0:
 63    if clearDirectories and os.path.exists(dataDir):
 64        shutil.rmtree(dataDir)
 65    if not os.path.exists(restartDir):
 66        os.makedirs(restartDir)
 67    if not os.path.exists(vizDir):
 68        os.makedirs(vizDir)
 69mpi.barrier()
 70
 71#-------------------------------------------------------------------------------
 72# Material properties.
 73#-------------------------------------------------------------------------------
 74units = MKS()
 75eos = GammaLawGas(gamma = gamma,
 76                  mu = mu,
 77                  constants = units)
 78
 79#-------------------------------------------------------------------------------
 80# Create our interpolation kernel
 81#-------------------------------------------------------------------------------
 82WT = TableKernel(WendlandC4Kernel())
 83output("WT")
 84
 85#-------------------------------------------------------------------------------
 86# Create a NodeList and associated Neighbor object.
 87#-------------------------------------------------------------------------------
 88nodes = makeSolidNodeList(name = "gamma-law gas points",
 89                          eos = eos, 
 90                          kernelExtent = WT.kernelExtent,
 91                          nPerh = nPerh)
 92
 93#-------------------------------------------------------------------------------
 94# Generate the initial node geometry (lay down the points).
 95#-------------------------------------------------------------------------------
 96generator = GenerateNodeDistribution2d(nRadial = nRadial,
 97                                       nTheta = 1,
 98                                       rho = rho0,
 99                                       distributionType = "constantDTheta",
100                                       rmin = rmin,
101                                       rmax = rmax,
102                                       xmin = (0.0, 0.0),
103                                       xmax = (rmax, rmax),
104                                       theta = 0.5*pi,
105                                       nNodePerh = nPerh)
106distributeNodes2d((nodes, generator))
107print("Point distribution across MPI ranks:")
108output("  mpi.reduce(nodes.numInternalNodes, mpi.MIN)")
109output("  mpi.reduce(nodes.numInternalNodes, mpi.MAX)")
110output("  mpi.reduce(nodes.numInternalNodes, mpi.SUM)")
111
112#-------------------------------------------------------------------------------
113# Set the node properties.
114# The above geometry initialization set the node positions, masses, mass
115# density, and H tensors.  We still need to initialize the energy for the Sedov
116# blast energy.
117# In this case we're going to simply pick the innermost ring of points closest
118# to the origin and dump the initial spike evenly across them.
119#-------------------------------------------------------------------------------
120pos = nodes.positions()
121mass = nodes.mass()
122eps = nodes.specificThermalEnergy()
123
124# Find the points closest to the origin.
125dr = rmax/nRadial
126spikePoints = [i for i in range(nodes.numInternalNodes) if pos[i].magnitude() < 0.75*dr]
127print("Selected a total of {} points for energy deposition.".format(mpi.allreduce(len(spikePoints))))
128
129# Distribute the energy evenly in energy/mass.
130Mspike = mpi.allreduce(sum([mass[i] for i in spikePoints]), mpi.SUM)
131eps0 = Espike/Mspike
132for i in spikePoints:
133    eps[i] = eps0
134
135# Verify we actually initialized the expected energy.
136# We are using multiplication of Spheral Fields here, and the fact that Field.sumElements
137# sums across all MPI domains by default
138Esum = (mass*eps).sumElements()
139print("Initialized a total energy of {}".format(Esum))
140assert fuzzyEqual(Esum, Espike)
141
142#-------------------------------------------------------------------------------
143# Construct a DataBase to hold our node list
144#-------------------------------------------------------------------------------
145db = DataBase()
146db.appendNodeList(nodes)
147output("db")
148output("  db.numNodeLists")
149output("  db.numFluidNodeLists")
150
151#-------------------------------------------------------------------------------
152# Construct the hydro physics object.
153#-------------------------------------------------------------------------------
154hydro = FSISPH(dataBase = db,
155               W = WT)
156packages = [hydro]
157
158output("hydro")
159output("  hydro.cfl")
160output("  hydro.compatibleEnergyEvolution")
161output("  hydro.densityUpdate")
162output("  hydro.HEvolution")
163
164#-------------------------------------------------------------------------------
165# Create boundary conditions.
166#-------------------------------------------------------------------------------
167xPlane0 = Plane(Vector(0.0, 0.0), Vector(1.0, 0.0))
168yPlane0 = Plane(Vector(0.0, 0.0), Vector(0.0, 1.0))
169xbc0 = ReflectingBoundary(xPlane0)
170ybc0 = ReflectingBoundary(yPlane0)
171for p in packages:
172    p.appendBoundary(xbc0)
173    p.appendBoundary(ybc0)
174
175#-------------------------------------------------------------------------------
176# Construct a time integrator, and add the one physics package.
177#-------------------------------------------------------------------------------
178integrator = CheapSynchronousRK2Integrator(db, packages)
179integrator.lastDt = dt
180integrator.allowDtCheck = True
181output("integrator")
182output("  integrator.havePhysicsPackage(hydro)")
183output("  integrator.dtGrowth")
184output("  integrator.lastDt")
185output("  integrator.dtMin")
186output("  integrator.dtMax")
187output("  integrator.allowDtCheck")
188
189#-------------------------------------------------------------------------------
190# Build the controller.
191#-------------------------------------------------------------------------------
192control = SpheralController(integrator = integrator,
193                            kernel = WT,
194                            restartStep = restartStep,
195                            restartBaseName = restartBaseName,
196                            restoreCycle = restoreCycle,
197                            vizBaseName = "Sedov-cylindrical-{}".format(nRadial),
198                            vizDir = vizDir,
199                            vizStep = vizCycle,
200                            vizTime = vizTime,
201                            SPH = True)
202output("control")
203
204#-------------------------------------------------------------------------------
205# Finally run the problem and plot the results.
206#-------------------------------------------------------------------------------
207control.advance(goalTime)
208
209#-------------------------------------------------------------------------------
210# Plot the final state if desired.
211#-------------------------------------------------------------------------------
212if graphics:
213    rhoPlot, velPlot, epsPlot, PPlot, HPlot = plotRadialState(db)

Now let’s go through each section of this script in some detail.