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.