[from Anthony] progress cont...

Just a quick comment on the asymptotic behaviour of D...what I thought for k3_real (red line) appeared as it looked was that we need to take care of the sign when doing a square-root...so one root for k=+sqrt(-D) and another for k=-sqrt(-D) which was why k3_real looked like a reflection of k1_real for r>>1. I agree...we need to do more analysis to explain k2 and also k3_real for r<<1.

****

So I tried integrating starting from the left with the following initial settings...

m=2 (yes...m was equal to 2)

dr = 0.003 (small enough to see if there could be kinks or interesting lumps)

r_initial = 0.01

eta(r=r_initial)= 0

v_eta(r=r_initial)=1

The coefficients B(r) and C(r) should be correct as I did the asymptotic predictions analytically and confirmed it numerically. Oh, from what I found, for r>>1, C(r) is actually increasing (rather than decreasing from our other parameters) because of the last term D/c^2, which is directly proportional to r (assuming c ~ v_kep). Nevertheless, the only part I think that has the greatest contribution to the solution of eta is I(r). Doing the asymptotic predictions for I(r) was a problem for me because I was not sure what assumptions I can do to simplify the term so I had to do it numerically.

Before I show the results, I'm going to define some of the symbols I used...

tau_1 ~ perturbed optical depth

c ~ sound speed

v_tau ~ dtau/dr (for unperturbed optical depth)

f ~ unperturbed radiative force

WW ~ shorthand for m*(Omega-Omega_p)

sigma ~ surface density

a ~ 2nd derivative of eta

v ~ 1st derivative of eta

From the list: sigma, c, f, WW, and v_tau have already been calculated prior to this integration.

So with the initial settings I adapt the leapfrog method for the integration...

tau_1[j+1] = tau_1[j]+(1/(c[j]*c[j]))*v_tau[j]*eta[j]*dr;

		

I[j+1]= I[j]+f[j]*(((log(r[j+1]*sigma[j+1]*f[j+1]/D[j+1])-log(r[j]*sigma[j]*f[j]/D[j]))/dr)+

				   ((2*m*W[j])/(r[j]*WW[j])))*tau_1[j];

		

a[j]= - B[j]*v[j] - C[j]*eta[j]- I[j];

		

v[j+1] = v[j] + a[j]*dr;

eta[j+1] = eta[j] + v[j+1]*dr;

So again, the only part I'm suspicious is tau_1 which includes I(r). But here's what I got from this....

I stopped the integration at R=5 because I thought it was meaningless to continue on. This seems to act very similar to the wave for k_3 (look back at my previous post). Now, I thought it was weird that 'nothing' was happening between r=0 to r=1 so I restricted r to this range and I get...

I noticed something was happening between r=0 and r=0.2 so I restricted r further...

Looking at the numerical results, it seems that C(r) was the dominating factor (which I think was the one making these responses) and after r=0.4, I(r) became the dominate one very quickly. Moreover, this also resembles the waves for k_1 and k_2. I also plotted eta between 0.4 to r=0.8-ish...

So if you piece together the last 3 pics you would get the first pic I posted.

That's all I have so far in terms of integrating. If you need to see more pics, let me know and I'll post them.

***

Pawel, I believe you reminded me that this second ODE we been dealing with is actually an integro-differential equation. So the only method that rings to my head to solve this kind is using Laplace transform. With the initial conditions we defined (eta(r_initial)=0; v_eta(r_initial)=1) it would be convenient to use this method where we transform and solve for L(eta) then inverse laplace transform it to get an analytic solution. The first couple of terms are easy to transform but for I(r) I was not sure whether we should treat (1/c^2)(dtau/dr) as a constant and pull them out of the integral. If I did that, it would be easier to solve. So now I'm at simplifying and solving for L(eta) and so far I have this fraction with a 3rd degree polynomial at the denominator. Now I just need to find some way to simplify it and inverse transform it. But before I go further, do you think this method is actually suitable for our case???

Anthony

[From Jeffrey] cuda-hydro

So the code works now...



The resolution for this simulation is 732x732 ( in case you wonder why such strange numbers: 732=3*(256-12) ). Approximately 10000 time steps were taken and it took 2 hours on k3. Total duration in simulation time is approximately 2 crossing time.

This is definitely not the fastest it can get. Computations are currently only done on one block that has a size of 256, so the simulation box is 3x3=9 times too large for the block. I should be able to utilize 9 blocks in parallel to do this simulation, which would in principle give me another factor of 9 in speed. The only reason why this is not done yet is that the memory sharing between blocks is not trivial and I didn't want to involve too much complication before I get it to work. Now that I understand cuda a little more, I already have a plan for the improvement and it shouldn't be too difficult.

[from Anthony] Progress report on modal analysis

First of all, I'm now back to Toronto from Santa Cruz and I have to say ISIMA was excellent. I definitely learned that there's a lot of more complicated problems in astro modeling where it involves lots of hydrodynamics (MHD) and there was one group from UC San Diego working on fusion (plasma physics). So, after meeting a bunch of graduate students from all over the world there's definitely a lot more we could learn in astro modeling.

Here's the site for what's currently happening in the program...

http://isima.ucsc.edu/current.html

***

So while I was at Santa Cruz, I managed to get some work done on my modal analysis with Pawel. Having to solve for the roots of k (wave number) via WKB approximation along with Cardano's method on the third order polynomial; so there are three roots to this polynomial.

This is the first root of k with both real and imaginary parts. Notice the the three lumps on the graph. Those are the Inner Lindblad, corotation, and Outer Lindblad resonances respectively.

This is the second root for k.

This is the third root for k.

By integrating for eta (look back on my previous entry on what eta is...) using the WKB solution exp(ikr) or more specifically exp(i*integral(k_{i} dr)) we get the following three eta's...

This is the first eta using k_1 with both real and imaginary parts. There seems to be a large response between 0.1 to 0.5 and then some response between r=1 to r=1.7-ish.

This is the second eta using k_2. It seems nothing happens after r=0.5.

This is the strangest one. It seems eta acts violently (look at the number on the eta axis) as r grows which I can't seem wrap my head around. Looking back at the roots of k, k3 sort of acts opposite to k1, so there's a sign reversal on the k's which implies that the graph is flipped horizontally but I don't quite understand why there's such large values for eta_3.

[from James] current status of dust-disk sim

Currently there are a couple of versions of simulation based on Josh's original code. On Cudak4, they're somewhere inside the student/.../C/src/ directories (I don't recall the exact path but it's the same folder in which all the other CUDA programs are kept). There are two versions in two separate folders called "tau" and "taunot", respectively.

For some reason, the program which ran just fine on cudak2 and cudak3 Sunday night wouldn't run on cudak4 or cudak5; the first two cudaMemCpy calls (from host to device) worked fine, but all subsequent cudaMemCpy calls, regardless of direction, threw a "wrong direction" error. Equally strange, commenting out everything tau-related made it work again.

The "taunot" folder contains the working simulation in which everything relating to the tau grid has been commented out (which is why, if you re-compile it, nvcc complains several times about nested comments). Also in this version, just for the fun of it we included a force term in the theta'' calculation to simulate a non-rotating bar. As mentioned previously, you can also ssh into cudak4 and run this sim remotely (although when I tried it from home the display was hideously slow, presumably because I have a slow internet connection and a 5-year-old computer).

The "tau" folder contains the non-working sim-in-progress in which the tau-grid is being re-implemented. After much discussion and debate (even some argument) on Tuesday we decided to re-implement tau using the following scheme:
- use two grids: an integer array "dtau_set" storing the number of particles in each grid cell; and a float array "tau_set" storing the value of tau in each cell (this will of course take up more memory but for now it seems that speed is more important)
- when the disk is initially generated, find out which cell each particle belongs to and increment that cell in the integer dtau_set
- start the simulation loop
- kernel 1: send one thread along each "ray" (solid angle) from r_min to r_max, calculating tau for each cell in the float array tau_set based on the data stored in the dtau_set integer array
- kernel 2: integration, one thread per particle: first figure out which cell the particle is in, retrieve the appropriate value of tau, and calculate beta; then calculate the acclerations and update the new positions and velocities; now re-calculate the particle's cell location, and if it has changed cells, decrement the old dtau_set location and increment the new dtau_set location

That's the plan, still in progress. Feel free to work on the code. Clearly this algorithm still isn't super-efficient: there are several calls to global memory, and some values get re-calculated more than once even though they haven't changed (most notably the particle's cell index). However we have managed to limit ourselves to two kernels, and each thread is independent with no need to invoke syncthreads(). The presentation date has come and gone (with no presentation :( ) but for my part I'd still like to see this thing up and running.