PHYD57 Advanced Comp Methods in Physics - A4 Projects in 2024

Introduction

Our course gives you an opportunity to delve into a subject of your choice (to some extent!).

You will study the subject, create algorithms and plan the visualization of results, implement them to accomplish some numerical task such as a scientific simulation, debug the code and run it either on your on our linux server art-1, under your subdirectory (there will also be 4 subdir's dedicated to 4 approved projects, for file sharing & co-developing programs).

Form a group

The assignment involves creating a writeup and giving a presentation (cf. below). There will be 4 groups. Forming groups (max. 4 persons ea.) is up to you; discuss it during breaks between the lectures. Choose group/topic by the time of midterm. . If you don't, a topic will be assigned to you after misterm. Send email to me (P.A) as soon as you know two ranked, preferred A4 topics. For instance, if 5 students choose some topic and only 1 some other topic groups will have to be rearranged. It is advantageous to talk to fellow students a.s.a.p. See list of emails hopefully protected from search engine bots & general public [I knew you'd click :-), but the page you see doesn't have sufficent protections; I linked the emails at the top of our../XX auxiliary resources page]. Send me the proposed list of names (can be less than 4) with the topic number. I will try to reserve the project for you on a first-come-first-served basis. Some ideas are listed below. The final, reduced list should be available by mid-October.

Marked components of the project

Show your initiative, curiosity and ambition, do interesting computations and tell us about them, even if they deviate from the suggested descriptions below. All group members get the same mark for A4, don't tolerate freeloaders.

1. Every group member works on and presents part of a common PDF or Powerpoint presentation. All group member should present part of the talk - everybody gets to practice that. Each group has about 25 min to present the work, followed by usually just a brief 5 min discussion (questions from audience) -- unless the topic raises unusual amount of interest or controversy. It's easy to go over the allotted time; try not to, by practicing and timing the talk in advance. Points will be subtracted if a group badly misjudges in either direction the time taken by the presentation.

2. Prepare 8 to 12 page PDF writeup (incl. figures), due 1 day before your presentation, so that we can read it and think of interesting questions that we ask you (everybody - try to prepare a question for the other projects!). Writeup may be partially overlapping with but is NOT the same as pt. 1 above. Writeup is like an essay, popular article or long blog entry. It includes: names and student numbers below the title, then some background on the problem, methods others and you have used, your code's details, tests showing that it actually works against all odds & returns meaningful results (provide simple test cases), followed by your results (estimate accuracy!). Writeup ends with a section giving analysis and then draws conclusions about the studied system. Unlike the slideshow, writeup needs to include the full list of references to cited work, including internet pages, if any. Slideshow - not necessarily all, only the most important ones.

Topics (# to be reduced to 4 or 5 by midterm):

Detailed description of most projects

1. Parallel simulation of planets named after a telescope, named after a beer named after monks

See L11 for background on chaos in dynamical systems, and the concept of Liapunov exponent (or inverse of it, called Liapunov time). See previous lectures on 4th order symplectic integrator that you are required to use. See all CUDA materials on how to write the CUDA kernels that concurrently integrate $K \gg 1$ slightly different systems (differing in initial positions and/or velocities of planets). The calculation should not be serial, i.e. multiple versions of a planetary system described below should be simulated concurrently on CPU or GPU. [The choice is yours, and may warrant a little experimentation first. Namely, it's not clear if a better overall throughput measured in number of simulated system-years per second of wall clock time, is achieved on CPU or GPU, because both calculations need to be done in double floating point precision & consumer GPUs on art-cluster are not as fast in double as in single precision]

Data for the fiducial (basic) setup are largely specified in the entry of the catalog http://exoplanet.eu/catalog (enter Trappist-1 into search box). Trappists are Belgian monks making a famous dark beer (consult LCBO, if interested). They did not observe exoplanets. Graduate students doing observations apparently learned a lot about the Trappist's products during their studies. The system has 7 planets, all of them Earth-class and a few in the habitable zone! Find out more about it on the net.
Perturbations to the initial conditions away from the given set of parameters will be small and pseudo-random (if they were *really* random, you would not be able to repeat the same numerical integration, which you need to do to debug and verify your code).

If you haven't got the knowledge of Kepler problem (2 Body problem) from ASTC25 course to set up initial conditions of the simulation, read about Kepler equation and its iterative solution on the web and ask me questions diring tutorials. Assume that at the start of the calculation, the phase of the motion (fraction of full orbital period the body did from the last pericenter), orientation of the orbit, and hence also the true anomaly $\theta$ (azimuthal angle in polar coordinates; it is figuring in the following, parametric, equation of ellipse: $r(\theta) = a(1-e^2)/(1+e \cos \theta)$) are all random numbers with uniform distributions (e.g., $0 \lt \theta \lt 2 \pi$ rad). Both statements cannot be simultaneously strictly true, since the motion on ellipse is not uniform in time. But since the Trappist planets are on fairly circular orbits, the assumption is justified in practice. Analytically, i.e. with pen and paper, from the equation of ellipse cited above derive the instantaneous radial and transverse components of velocity vector at any theta, and recalculate to Cartesian coordinates if you use that system for calculations. A hint is that angular momentum per unit mass is a constant in the elliptic motion: $L = r^2 \,d\theta/dt = r v_\theta = const. = \sqrt{GMa(1-e^2)}$. Therefore $v_\theta$ follows very easily from the initial r and the parameters of the ellipse, and $v_r$, which contains derivative $d\theta/dt$, is also easy to compute.

There are actually two ways to quantify the stability, one by watching for close encounters and/or escapes and finding the time at which the change occurs, the other by integrating the shadow systems with miniscule initial deviation of the positions and velocities, to find the Liapunov exponent. I would like you to use the first one and the 2nd as additional indicator if at all. We want to label a system "unstable" if at any point in time (you can check every time step or every now and then, not as often as every time step) planets approach to within 2 Roche lobe radii (see L12) of any of the approaching planets, or if a planet is located further than 2 times the initial distance between the star and outermost planet, or closer to the star than the initial radius of the innermost planet. Ideally, the calculation would then stop and a new variant of the system started. But it's up to you how to optimize & parallelize. Describe everything in the writeup. Also, in any dynamical calculation describe how you chose the optimum constant time step dt, and how you tested the accuracy of your simulation.

2. N-body galaxy with a migrating Black Hole

This is essentially an N-body simulation with N on the order of 100 million, in which N stars in a target/host galaxy and in a smaller dwarf galaxy (that has just collided with the host galaxy) are feeling the attraction of all the other stars. Direct computation of all the interactions is too slow. All gravity forces needs to be evaluated on a 3d grid by FFT method already outlined in L10 (more details in L12). Zero-padding will be used to avoid aliasing; the size of full density array will be 512x512x512 or at least 512x512x256 (flat host galaxy and nearly coplanar initial orbit of the dwarf one allow a smaller resolution in "z" direction). Each of the two galaxies hosts a relatively massive particle in the center, representing an SMBH (supermassive black hole).

The objective is to write and test the FFT gravity solver first, then test the dynamic integrator (2nd order symplectic, leapfrog), show that both work accurately enough (time step for dynamics integrator must be optimized). Testing involves checking the conservation of those physical quantities that should be constant. Finally, you'll run the code for an extended period of the simulated time, and see what happens with the host galaxy, dwarf galaxy, and their SMBHs. There is no limit how much fun you can have doing the visualization. The minimum is to produce a series of snapshots in 2-D projection(s); plan maximum would be to do animations/videos.

3. SPH project

See lectures 10 and 11 on SPH, both nearest neighbor finding using linked lists, and the SPH algorithm to do gas dynamics.
The best description for the purpose of implementation that I've found is this article in Annual Review of Astronomy and Astrophysics:
Smoothed Particle Hydrodynamics (1992).
In addition, here you have two other descriptions with more details, derivations and comments on implementation:
Monaghan (2005).     Monaghan (1997).
Description from Lagrangian point of view (optional material! I frankly can't remember which book it is from) Chapter 3 (2010).
You study aresponse of to a gravitational perturbation of a disk of gas. In principle, you can either study any disk such as a disk of gas in a spiral galaxy, perturbed by an imposed, known, low-amplitude spiral disturbance of gravity from a stellar disk (which you do not need to model, though in real research you would have to - in a way explained in the previous project!). I suggest that you study a protoplanetary gas disk disturbed by the gravity of an embedded protoplanet. You should program the planet to travel on a circle at a speed described by Kepler's laws, but introduce it into calculation gradually, leaving time for the disk to relax both (i) the initial 'ringing' due to imperfectly calculated initial conditions of force balance, as well as (ii) the unavoidable oscillations purely due to rapid insertion of a perturber into the disk. Therefore, ramp the mass of a planet linearly up to the full mass (study 2 mass ratios = m/M = 0.001, 0.01) in time t=0...10 orbital periods. The planet should be in the middle of the disk at radius 1 (if you prefer to use correct physical units instead of rescaled non-dimensional coordinates, then 1 would mean 1 AU = 150e6 km). The disk should extend from 0.45 to 2 (AU). On demand I'll describe the disk setup in more detail. (The disk should be 3-D, and have z/r = 0.05 throughout, where z = vertical scale height.) Mass of the disk should not play a role. Disk should be vertically isothermal, its soundspeed should only vary with the cylindrical radial coordinate r. What maximum density contrast in the disk do you obtain after 10, 20, and 50 orbital periods of the planet?

4. Pilot's glory. Mie theory calculations.

Pilot's glory (a.k.a. Brocken Spectre; cf. wiki article on glories) is a colorful halo surrounding the shadow of an object (observer) cast on the cloud or fog consisting of suspended water particles. It results from a near-180 degrees backward scattering of sunlight off a collection of spherical water droplets in the air. (That's why it is seen surrounding a shadow).
Picture 1    picture 2
It can be predicted. And vice versa: from observation, in principle, we should be able to find out the distribution of droplet sizes corresponding to each picture. You can find more pictures here, as well as here , and choose one to explain.

In the phyd57/mie subdirectory you will find some working fortran subroutines such as mie.f90, predicting the spatial distribution of light scattered off a spherical water drop, for a given optical refraction index (or alternatively dielectric constants) of water, and ratio of the particle size to the wavelength of light. [There are in fact 2 separate inputs, radius of a particle and the wavelength of light, since the refractive index of water at a given wavelength can be computed by a separate little function or procedure.] Using that subroutine or translating it into a C(++) procedure, you will assemble a theoretical color image of pilot's glory using 200 separate wavelengths of visible light, and 100 different sizes of water droplets. (You can use Planck function with temperature 3700 K to model the spectrum of a slightly reddened sunlight.) The distribution of water drop sizes will be Gaussian, characterized by a mean size and half-width of distribution (total width will be +-3 half-widths). Fiducial mean size is lambda = 10 micrometers, and half-width of the distribution function is delta = 3 micrometers. Ideally, you will be able to adjust the two parameters of the size distribution until you find a reasonably good match to the chosen picture. At present, nobody knows what size distribution the cloud shown in the picture had. You may be the first person to know it.

If you are interested in Mie theory, although this is not strictly required to master for this project, here you have a book about it. There is a fortran code in the appendix, though my old code is a bit better :-)

5. Machine Learning from very noisy data

Design and code your group's machine learning program. Make it solve some tough problem in econophysics and/or algorithmic trading. Using historic data, teach the machine to buy and sell stuff in a time-variable market, at least to the extent given by the short time-frame of this projects. Extra points will be awarded for making your robot crack the secret of the stock market, and paying off your and your fellow students' loans. With skills in Python and a little patience, high-resolution Wall Street historic stock data (1 minute resolution) can be downloaded in chunks for free. As you learn about machine learning, you will likely realize why lots of data are needed to teach robots non-trivial things.

6. Structural dynamics: build and dynamically load a virtual Golden Gate bridge or Citicorp Center skyscraper

There are many stories and books about Golden Gate. And about structural engineering. Make a digital twin of a bridge that responds to dynamically moving loads on the roadway, study the deformation and vibration modes. Did you know that during one anniversary the traffic was closed and so many people gathered on the bridge that it visibly deformed, and the police had to dissolve the celebration in a hurry? How many people are too many?

Did you know that a Citicorp Center skyscraper in NYC was realized independently by two students (doing something like your project) to have a possible, fatal flaw? It could be catastrophically unsafe in a wind exceeding just 110 km/h. An approaching Atlantic tropical storm carried such winds, but passed Manhattan safely. In great secrecy, Citibank tower was reinforced while in use. Tell us the story, but mainly build a model of a skyscraper and load it with aerodynamic stresses due to wind gusts. Analyze results, test how realistic your model parameters are. How high a wind collapses your skyscraper?

7. What happens to an airliner after it loses a part of the wing: panel simulation of ideal aerodynamics of airfoil

Study what happend on 10.04.10 in Smolensk, Russia, to a Polish presidential airplane, by physical calculation. Learn a panel method of interaction of an air flow with an airfoil, which is different from SPH and grid methods of CFD (Computational Fluid Dynamics). It uses Kelvin's circularion theorem, Amper's law, Linear algebra, and represents the 3-D flow in infinite space by solving integro-differential equations. Start by reading great 100-year old books that are still used today in preliminary design of aircraft! Using the calculated aerodynamical forces, predict the trajectory and orientation of the airliner after collision with a large tree and compare with accident reports. Confirm or debunk conspiracy theories surrounding this catastrophe.
See sources on our auxiliary /xx page. You should get acquainted with the great 100-year old research by a Bavarian physicist Ludwig Prandtl, who created whole school od aerodynamics, both theoretical and experimental. His students included Theodore von Karman, who became the father of American theoretical aerodynamics & helped U.S. develop first supersonic jets, and Richard von Mises, who achieved many outstanding results in elasticity theory of solids, boundary layer theory and also airfoil theory. Especially von Mises needs to be mentioned in the context of bad political ideas Prandtl upheld during the Nazi regime in Germany in 1933-1945, which Prandl generally liked and praised (Mises was a Jew from Lemberg = Lwów = L'vov, now Lviv).

Page through Prandtl's books
Prandtl - Funadamentals of Hydro- and Aeromechanics
Prandtl - Application of modern hydrodynamics to aeronautics
The more modern re-telling of classical works by Prandtl is the subject of
Oertel et al. - Prandtl's "Essential Fluid Mechanics" (2004)
The part you are particularly interested in is Lanchester-Prandtl lifting line theory of airfoils.

Your task is to first create a lifting-line theory of a full (test case) and a damaged wing of TU-154 airliner (for which you'll only need the top view of the aircraft and the distance from its axis where a birch tree tore off about 1/3 of the left wingspan, and some infor about the AOA = angle of attack). The theory is formulated as in integral equation, discretized it becomes an NxN linear equations problem, which can be solved with Python. It yields the dirtribution of stationary lift force along the wing of any shape and asymmetry, so you can derive a torque on a damaged aircraft, needed to predict the precise history of the roll to the left that develops, wheter or not the aircraft was controllable by pilots at this stage and such.

Next, if you have enough time, you will generalize the method to full panel method in which vortex lines are layed out not in the form of lines like in Prandtl's scheme, but around N polygonal panels covering the upper and lower surface of the wing. I haven't gone so far in my calculations, we'll see if the panel theory gives similar results.

(Most of the above links are dead. There are live links on our /xx page and on the hidden pages described there)

8. Other topics

You consult with your friends and propose a topic. It will be pre-approved by the lecturer to assure appropriate complexity level and relevance to physical science. Say, a group may solve a wooden puzzle with 25 wooden pieces by exhaustive search of the tree of possibilities, and be rejected based on too little overlap with the subject of course. Or try writing a neural net program that controls some physical process or experiment (of course, that would be a simulation not a real-world experiment). Neural nets are a growing area of interest in physical sciences, so the project might be approved.

Quantum computing is an interesting topic, but how will you do computations using or related to it & what is the connection to numerical analysis of the real world?

One such self-defined project has been approved: Searching for exoplanets in HARPS spectrograph/telescope archive data. Contact Daniel Harrington.