The SFT Code: Finalised!

I finally have good versions of the code, which have passed every test I know!

One version is very slow, but very accurate. It runs with a time step of 0.9 seconds only, which gives us a runtime of 325 days for a simulation of 11 years i.e. one sunspot cycle. This code can be used for observing small scale features which do not exist for long.

The faster version of the code runs with a time step of 250 seconds, which gives us a runtime of 24 hours per sunspot cycle. In this code, the grid points within 4.32 degrees of latitude from the pole have 50 points in the toroidal direction, whereas the rest of the latitudes have 1000 points in the toroidal direction.

The evolution of the radial field is second order accurate in both forms of the code, and there is a slight flux imbalance due to numerical round off, which is less than 0.1 % of the total flux.

I ran a simulation with the modeled sunspot input for one solar cycle, and the output is shown in the video here.

A Working Code with a New Grid

The Grid used for calculations in all previous codes was:

theta= 0 to pi in 501 steps, with a step-size of h=pi/500.
i.e. theta=linspace(0,pi,501)

This grid includes both the poles (theta=0 and pi) and the equator as grid points.

The corresponding area grid took into account the area surrounding these points. i.e. the area between theta-h/2 and theta+h/2.  So, the area corresponding to the point just after the pole would lie between the colatitudes of h/2 and 3h/2, and would look like a trapezium. Whereas, the area corresponding to the polar grid point would lie between colatitudes 0 and h/2, and would look like 2 triangles with a common vertex, and a common axis of symmetry. And the area-grid of the points on the equator lied equally in both the hemispheres.
The main problems we faced due to this grid, were that

1. there were terms in our master equation, that involved 1/sin(theta). These terms blew up at the poles.
2. The problem of flux transport over and across the poles (singularity).

Solution: Shift the grid by h/2.

The grid was shifted by h/2, such that it did not include the poles or the equator, or any integral multiple of pi/2. The new grid is:

theta= h/2 to pi-h/2  in 500 steps, with a step-size of h=pi/500.
i.e. theta=linspace(h/2,pi-h/2,500)

The corresponding area grid, includes triangular area for the points next to the pole, but now there is only one triangle on one side of the pole between colatitudes 0 and h, so we don't have to worry about flux crossing the poles. (Because the transport of flux in the area grid occurs between 2 areas with common side. But at the pole, the contact between 2 areas in the theta direction is through a single point, which gives us no flux transport across the pole. This is also theoretically verified.).

And the 1/sin(theta) terms in equations do not blow up because there is no grid point exactly at the poles.

So, our purpose of solving the issues is served with a shifted grid.

After doing this, another problem which is right at the heart of all computational tasks remains!

Computing time:
Initially, upgrading the code with the new grid, a trial run was set, which had an estimated runtime of 2 years!

Then after some optimization and modification, I was able to reduce it to 1.75 years. Still not enough to even try a trial run this year!

So I looked at which task requires the least time step. It was the diffusion in the toroidal direction. It requires a ridiculously small time step of 0.5 seconds! This is because the grid points near the pole are very close.

So, I did the following:
Made 10 groups of 100 grid points in toroidal direction (total 1000 points on every latitude) for the first 10 colatitudes from the poles.

Now, I calculated the maximum time stepping I could use for the rest of the grid except the first 10 latitudes. It came out to be almost 400 seconds. (Much much better! 800 times!)

So now, the code treats the first 10 latitudes explicitly, and evolves their magnetic field 800 times with a time step of 0.5 sec, and then taking that field as the boundary condition, it evolves the field at rest of the points once with a time step of 400 sec.

This reduced the runtime to 16 hours per solar cycle. :)

The very first run with new grid:
download the video or watch in HD here. 

Cons:

The flux transport is the same, as calculated from previous simulations with approximate boundary conditions. But the peak magnetic field is lower than expected.(Same flux, but low magnetic field) May be a result of unnecessary diffusion due to first order accurate scheme for meridional flow.