Set #4 of assignments in PHYD57, 2022.
=======================================
15 % of course mark (more than previous sets)
Problem 4.1 [30 pts] - Irradiation Instability simulation
__________________________________________________________
Linked to the programming page of our course (../progD57) is a program
tau-dis-cpu.f90 (also in dislin directory at phyd57@art-2.utsc.utoronto.ca)
With the help of OpenMP on CPU, this program simulates a disk of 1 M particles
around a star that exerts both gravity and radiation pressure with coefficient
beta = 0.2. The physical problem was discussed in tutorial T7 on day 9.
It consists of two parts: one is the Newtonian dynamics of N particles,
each of which orbits a point mass (particles are generated in a wide ring).
The second part that can be tested separately is connected with radiation
pressure on particles. Acceleration equal to gravitational acceleration times
1-beta*exp(-tau)
where tau = optical thickness = dimensionless number specifying what fraction
of the star (as seen from the particle's position) is covered by particles
in a in front of the given particle. Even if tau > 1, i.e.
combined area of particles is larger than the cross section of a channel
traversed by a stream of photons, there is always some transmitted light,
because the obscuring particles partially overlap. See wikipedia on optical
thickness and the exp(-tau) transmission law.
OBJECTIVE
The current version of the program partially works but not fully correctly.
Develop or rewrite & debug it in your favorite language (not Python).
DISLIN graphics commands are virtually the same in both C and F90 languages.
DO NOT WORK ON THE SOURCE FILE PROVIDED. Copy the program to your own
subdirectory and rename to include your first name or 3 digits of
student number, before changing.
PHYSICAL SETUP
Gravity (in code units. where central source of gravity gas GM=1) is
f_grav = -1/r^2
and the radiation pressure acceleration is
f_rad = +beta*exp(-tau(r,phi))/r^2
where beta = 0.2 is a constant (relevant is there the full starlight
falls on particle) and exp(-tau) is the attenuation factor that
is computable from the density distribution of particles in front
of a given radius (i.e. on the line of sight between the star and the
observation point).
tau(r,phi) = optical thickness = integral from 0 to r { const/r^2 dr}
const. = c * Dirac delta.
What do we mean by this? What is meant here is this. Each particle, by
assumption, has the same cross sectional area. Consider a conical disk sector
around the azimuthal angle phi=const., which also grows linearly in height
(physical disk is rather flat, but assuming strict 2-D structure is not
correct). The cross-sectional area of the sector at radius r grows in
proportion to r^2. Each particle contributes to optical thickness an amount
which is just the ratio of the its cross-section to the sector cross-section:
~1/r^2, or more precisely: const./r^2.
The particle is localized in radius and formally looks like Dirac delta
function, so const. = c * Dirac delta(r-R)
(peak of Dirac delta is at particle's radius R, c is a positive constant).
Substitution turns integration from a continuous integral to a sum over all
such particles in a given phi-sector, which are closer than r to the star.
The value of constant c is computed by the program in the first
call to the subroutine rad_pressure like so: c is set to 1.0 at first,
which incidentally allows to measure the very small collision rate of
updates to histogram array dtau during the multithreaded operation of
the program, and in all later calls c is kept unchanged. Its value is
fixed at whatever it needs to be to assure that the total optical
depth across the entire initial disk out to Rout=2.5 radius is equal
tau_max=5.
This makes the disk optically thick, i.e. non-transparent, and in turn causes
what we call IrRadiation Instability. IRI was discovered in gaseous disks
and described in 2015 in Astrophysical Journal by you lecturer and a former
UTSC student Jeffrey Fung (now asst. prof. at Clemson Univ.).
IRI also happens in purely particle disks (see animation in L9 recording).
Angular momentum of each particle is L=const., since there are no tangential
forces in the problem, all forces are radial.
We compute and store L values in slot 5 of the particle data array partic(6,:).
In fortran90, ":" is like asterisk *, it means all values of the second index
(particle number i=1...Np).
Transverse speed does not have to be found by integration, since it directly
follows from the definition of L: v_phi = L/r.
The radial equation of motion involves physical accelerations
[beta*exp(-tau)-1]/r^2 [radiation + gravity], and the centrifugal force
f_cent = v_phi^2/r = L^2/r^3.
There is no friction or collisions, so
d(v_r)/dt = [beta*exp(-tau)-1]/r^2 + L^2/r^3
v_phi = L/r
is a pair of equations (one differential, one algebraic) from which velocity
components follow, and from them one can easily find the positions (r,phi)
by integration over time:
dr/dt = v_r
d(phi)/dt = v_phi/r = L/r^2.
If tau=const, then these equations result in a closed, elliptic trajectory.
You can check whether the code produces such a trajectory by setting tau to
a constant and plotting a few of them.
WHAT TO IMPROVE
First, if you'd rather write your own code, do it. Sometimes re-writing
parts of a code from scratch is a good way to remove hard-to-find errors.
The Dislin graphics needs improvement. Right now, when you run
it on art-2, labels come out a bit funny. Links to documentation for
dislin are posted on our programming page.
Next and more importantly, there is something wrong with the
algorithm, because the expected instability does not seem to happen!
The disk remains axisymmetric.
You need to go into the code and insert temporary print instructions
in different sections of the code, in order to find which variables
have incorrect values, or if there is some mistake in the logic of
the program. For instance, there could be a situation where tau(r,phi)
is miscalculated, and for all phi-angles the r-dependence of tau is
not the expected function, which should rise from 0 to the value about 5 in
the outermost radius (index about 90% or 95% of the limit Nr in the code).
This would be a good hypothesis because IRI is caused by variable tau,
and the basic rotation of the disk seems working ok in the provided code.
Problem 4.2. [20 pts] - Theoretical image of beta Pictoris disk
________________________________________________________________
This problem is described here
https://planets.utsc.utoronto.ca/~pawel/bpic.html
and can be solved in one of: C(++),CUDA C, F90, CUDA Fortran.
Suggested graphics is DISLIN, but dumping file and producing a picture
in other graphics package or a Python script is allowed.
The essantial part of the code is a fast integration procedure.
This is needed, since a large number of theoretical images would have
to be created during an automatic optimization of the numerous
model parameters.
On request, more explanations will be given in the next tutorial.
Problem 4.3. [20 pts] - Linked list for 2-D close neighbor search
__________________________________________________________________
This problem can be solved in C(++)/F90/Julia (Python solutions will be
accepted but penalized with -25% of the mark).
Generate by a pseudo-random number generator described below,
N=30000 positions of points (particles) on a square of unit side.
Write a procedure that quickly finds the id's of all the neighboring
particles physically closer than h = 0.02 to particle number k (just
one particle).
In L9 and tutorial T9 a simple way to find all neighbors closed than h
was outlined.
Imagine a square grid with cell size h by h, and a catalog the
id numbers of each cell, using a linked list of length N (technically
implemented as a singly linked list). Having that liked list (integer
1-d array or vector of length N) and the grid array pointing to the
index of any of the points inside 'its' cell, which serves as the 1st
particle in a given cell, allows you to check the contents of the
nearest 9 cells. Not all particles will have physical 2-d distance < h.
That needs to be verified for each possible neighbor separately.
Generate N pseudo-random series of points on x axis this way.
Set a seed value to x = 0.123. Iterate inside a loop:
x = (12.64651+x*172.861); x = x - floor(x);
2N times to get (x,y) positions of N particles. Generate all x's before
all the y's, so that we all have the same positions.
The second instruction above leaves only the decimal fraction of a float
number. floor(x) rounds x down to a nearest integer. Such syntax should
work both in C(++) and Fortran90.
One needs to check the possible neighbors from only the 9 surrounding
cells, instead of all (1/0.02)^2 = 2500 cells, as would be true of a
simple-minded 'check-everybody' algorithm, trying (N-1) pairwise distances.
What is the computational complexity of your algorithm?
How many neighbors did you find for the point number k=100?
_____________________________________________________________________________
METHOD OF SUBMISSION:
Create a text or docx or pdf file with the output of your programs and
additional explanations, and submit to Quercus by the deadline.
In your subdirectory on art-2, create a subdirectory called A4 by $ mkdir A4.
Deposit working, well commented programs there.
Command line with best compilation instructions should be placed in the header.