Physics 581: Compute Anything –
Effective Scientific Computing#

Who: Graduate students or upper-level undergraduates from across the university with a strong preparation in mathematics (strong calculus and linear algebra with exposure to differential equations). Programming experience would be helpful, but we will provide support those with weaker backgrounds.
When: MW(F)2-4pm: Fridays are part of the iSciMath Coffee Hours and will be for working on problems, discussions, etc.

Imagine designing a rocket: You need to minimize weight, but maintain strength. Do you need to start a million-dollar program to model structures, or can you use computer modeling to help? What will get you to space fastest, and keep you in one piece? Compute Anything will lower the barrier for computing, allowing you to quickly explore problems as both a programmer and a scientist. We will show you that programming can fun, dynamic, and interactive, accelerating your research and your career.

If you have any questions, please contact Michael McNeil Forbes m.forbes@wsu.edu.

Programming Examples#

Here are a few examples of the types of problems we will solve in this course. Each of these examples is code-complete, using only NumPy and SciPy to compute the solution, and Matplotlib to plot (with a custom function FPS() to help make movies). By the end of the course, you should be able to write similar code within a matter of hours to study problems relevant to your research.

Molecular Dynamics#

Complex material properties can be simply modeled with molecular dynamics simulations that simply apply Newton’s laws to a collection of particles. Here we show a simple example of a projectile impacting a slab. See Molecular Dynamics for more details.

This idea can be easily generalized for your research: e.g. orbiting planets, galaxy formation, traffic flow, protein folding, and biomechanics.

Show code cell source

Hide code cell source

#%%time
import numpy as np
_EPS = np.finfo(float).eps

class MolecularDynamics:
    L1 = 250
    L2 = 40
    epsilon = 5
    sigma = 1
    m = 1
    N1xy = (40, 40)
    N2xy = (160, 40)
    v1 = 0 - 10j
    v2 = 0 + 0j
    r_cut_sigma = 2.5
    dt = 0.00005
    
    sep = 0+5j  # Initial separation (units of dx)

    # For colouring
    kinetic_energy_max = 5
    cmap = plt.colormaps['viridis']

    def __init__(self, **kw):
        for key in kw:
            if not hasattr(self, key):
                raise ValueError(f"Unknown {key=}")
            setattr(self, key, kw[key])
        self.init()
    
    def init(self):
        self.dx = 2**(1/6) * self.sigma   # Equilibrium spacing
        self.z0_1 = (-self.N1xy[0]/2 + (self.N1xy[1]-1)*1j +  self.sep) * self.dx
        self.z0_2 = (-self.N2xy[0]/2 + 0j) * self.dx
        
    def F_r(self, r2):
        """Return F/r where F is the Lennard-Jones force for separation r.
        
        I.e., F_ij = (F/r)*r
        
        where r = r_j - r_i.
        
        Arguments
        =========
        r2 : array-like
            Square of the separation distance.
        """
        # Prevent divide-by-zero errors
        r2_ = r2 + _EPS
        # s = (sigma/r)**6
        s = np.where(r2 > 0, (self.sigma**2 / r2_)**2, 0)
        s *= s**2
        f = 24 * self.epsilon * s / r2_ * (1 - 2*s)
        return f
        
    def get_initial_state(self):
        (X1, Y1), (X2, Y2) = [
            np.meshgrid(*[np.arange(N)* self.dx + x0
                          for (N, x0) in zip(Nxy, (z0.real, z0.imag))],
                        indexing='ij')
            for Nxy, z0 in [(self.N1xy, self.z0_1), 
                            (self.N2xy, self.z0_2)]
        ]
        Z1 = X1 + 1j*Y1
        Z2 = X2 + 1j*Y2
        V1 = 0*Z1 + self.v1
        V2 = 0*Z2 + self.v2
        Z, V = [np.concatenate(list(map(np.ravel, A)))
                for A in ((Z1, Z2), (V1, V2))]
        return (0, Z, V, 0*Z)
    
    def plot(self, state, ax):
        t, z, v, a = state
        kinetic_energy = self.m * (v.conj()*v).real / 2
        ax.scatter(z.real, z.imag, c=self.cmap(kinetic_energy/self.kinetic_energy_max))
        ax.set(title=f"{t=:.4f}", aspect=1)
        return ax
    
    def get_F(self, z, v):
        """Return the total force on particles at z."""
        r_ = z[None, :] - z[:, None]
        return (self.F_r(r2=abs(r_)**2) * r_).sum(axis=1)
        
    def step(self, state):
        """Return the acceleration on each particle.
        
        This uses the Velocity-Størmer-Verlet method.
        """
        dt = self.dt
        t0, z0, v0, a0 = state
        z1 = z0 + dt * v0 + dt**2/2 * a0
        a1 = self.get_F(z1, v0) / self.m
        v1 = v0 + dt / 2 * (a0 + a1)
        return (t0 + dt, z1, v1, a1)


md = MolecularDynamics(
    v1=-5j, N1xy=(4*2, 4*2), N2xy=(16*2, 4*2), sep=-10.3+2j,
    kinetic_energy_max=10, dt=0.00005*40)
s = md.get_initial_state()

fig, ax = plt.subplots(layout="constrained")
for frame in FPS(100, fig=fig, embed=True, filename="_static/MD.mp4"):
    for n in range(20):
        s = md.step(s)
    ax.cla()
    ax = md.plot(s, ax=ax)
    ax.set(xlim=(-12*2, 12*2), ylim=(-5, 10*2))
    ax.set_axis_off()
    fig.colorbar

Shallow-Water Equations#

Here is another complete example of one of the core applications that will be taught in the course: solving the Navier-Stokes equations in 1D for incompressible fluid flowing over a rough surface. In this example, we excite a fluid in a harmonic trap with a bump. As the fluid flows back and forth over the bump, viscosity dissipates energy and the sloshing slows down.

This idea can be easily generalized for your research: e.g. economic models, quantum mechanics, and biological response.

Details for those interested: The Shallow-Water Equations

The relevant equations follow from Newton’s law \(F=ma\), where the acceleration \(a = \ddot{X}\) of a particle at position \(x\) and time \(t\) is due to the external potential \(V(x)\), an internal energy \(gn(x,t)\) related to the hydrostatic pressure, and a viscous damping term with coefficient \(\nu\):

\[\begin{gather*} \ddot{X}(t) = \frac{F}{m} = -\frac{\nabla \bigl(V(x) + gn(x,t)\bigr)}{m} + \nu \nabla^2 \dot{X} \end{gather*}\]

This corresponds to an Eulerian hydrodynamic description of a compressible fluid with particle-density \(n(x, t)\) and velocity \(u(x, t)\) through the following partial-differential equations (PDEs)

\[\begin{align*} \pdiff{n(x,t)}{t} + \pdiff{n(x,t) u(x,t)}{x} &= 0\\ \pdiff{u(x,t)}{t} + u(x,t)\pdiff{u(x,t)}{x} &= -\pdiff{}{x}\frac{V(x) + g n(x,t)}{m} + \nu\pdiff[2]{u(x, t)}{x} \end{align*}\]

that can be solved by the method of lines as a system of ODEs integrating forward in time.

We will start with a simple implementation using finite-differences and Euler’s method to get an intuition for what is happening and to check our results:

\[\begin{gather*} \pdiff{f}{x} \approx \frac{f(x+\delta x) - f(x-\delta x)}{2\delta xh}, \qquad f(t+\delta t) \approx f(t) + \delta t \dot{f}, \end{gather*}\]

but then refine our solution using pseudo-spectral methods for derivatives and adaptive integration to get high-precision solutions. This latter solution is completed coded below to give you a flavour of the type of code that will be developed in the course.

There is some physics needed to describe details about this problem, like how the density \(n(x, t)\) relates to the positions of the particles \(X(t)\) in the first description, and how the internal energy \(gn(x, t)\) is related to the hydrostatic pressure. For details, see my Viscosity Notes 3 and Viscosity Notes 4.

The code here demonstrates the use of a pseudo-spectral Chebyshev basis to solve the shallow-water equations for an oscillating fluid. This code is complete, using only NumPy and SciPy libraries, except for plotting, which uses Matplotlib and a custom function FPS to facilitate making movies.

Not included here, but part of the class, will be supporting code for testing, additional documentation, and better organization to facilitate reuse.

For another example, see Model Fitting.

Learning Outcomes#

In addition to the Specific Skills listed in the Syllabus, by the end of the course, students should be able to:

Produce publication quality figures and animations (e.g., 2D and 3D plots, animations, tables, etc.) following best practices of visualization of quantitative data (i.e., as espoused by E. Tufte.)
Produce a manuscript meeting all formatting requirements for a paper in their field of focus, or for a thesis. (E.g., using a suitable LaTeX template and conforming to the submission requirements for mathematics, figures, citations, typesetting, etc.)
Work with data in a variety of formats for analysis, sharing, and reproducible science.
Use version control system to managing code or documents, and establish a portfolio of code on sites like GitLab or GitHub with testing via continuous integration (CI) and documentation.
Profile and optimize their code.
Use a code editor that supports regexp search-and-replace, linting, checking, etc. (E.g., Emacs, Vi, VS Code).
Numerically compute integrals and derivatives with quantified error estimates.
Use numerical linear algebra libraries to diagonalize and factorize matrices: QR, SVD, Cholesky.
Generate pseudo-random numbers with any desired distribution.
Use MCMC techniques to characterize and quantify Bayesian posteriors for fitting and model selection.
Use local optimization techniques (gradient descent, Nelder-Meade, BFGS, Broyden) with and without gradients.
Use techniques like Runge-Kutta to solve initial value problems.
Derive scaling results to validate these techniques.

Notes#

This is a live working document hosted on Read The Docs that will be used to collect and display additional information about the course, including:

Syllabus
Assignment Instructions
Resources, Readings, and References and various class notes. These should also be available through the navigation menu on the left (which might hidden if your display is not sufficiently wide).

These documents are built using JupyterBook (see JupyterBook Demonstration) and include all of the source code needed to generate the figure, plots etc. For example, to see how a figure was made, look in the preceding code cell. The complete source code for this documentation is available at wsu-courses/physics-581-physics-inspired-computation.

Funding Statement#

Some of the material presented here is based upon work supported by the National Science Foundation under Grant Number 2309322. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.

Physics 581: Compute Anything –
Effective Scientific Computing

Contents

Physics 581: Compute Anything –
Effective Scientific Computing#

FAQ#

Programming Examples#

Molecular Dynamics#

Shallow-Water Equations#

Learning Outcomes#

Notes#

Funding Statement#

Physics 581: Compute Anything – Effective Scientific Computing

Contents

Physics 581: Compute Anything – Effective Scientific Computing#

FAQ#

Programming Examples#

Molecular Dynamics#

Shallow-Water Equations#

Learning Outcomes#

Notes#

Funding Statement#

Physics 581: Compute Anything –
Effective Scientific Computing

Physics 581: Compute Anything –
Effective Scientific Computing#