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

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.

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

- there were terms in our master equation, that involved 1/sin(theta). These terms blew up at the poles.
- 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.