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.

One problem found in the simplest Advection term: Differential Rotation

The last term in the Master Equation,

is the advection term in the toroidal (phi) direction. This term is dealt with the Lax-Wendroff Scheme in the complex (and more accurate) code.

Simple first order accurate Euler scheme goes haywire after a few iterations, and does not conserve the peak of magnetic field (due to numerical diffusion).

But, the Lax-Wendroff scheme, which is second order accurate improves over the Euler Scheme, and preserves the peak of magnetic field, with negligible numerical diffusion.

A problem however, generally found with Lax-Wendroff schemes is, that some ripples can be seen behind the travelling wave. The total flux, and the magnetic field are still accurately conserved. But, due to the ripples, the unsigned flux suffers change with time; which is unphysical.

 The simulation of a meridian-like flux tube shows the formation of these ripples.

For a movie of the simulation, play the video named Flux tube1.mp4 after clicking here .

Problem with the Euler Code

While the simulation in last post had sunspots at 22.5 degrees latitude, and their flux was not carried over 80 degrees, the time step of 100 seconds gave us quite accurate flux conservation.

But, for a sunspot placed very near the pole, we need a time step of 0.1 second for a stable numerical evolution. When run with a time step of 0.1 second, the code will approximately run for about 3 years on my laptop to complete the simulation of one solar cycle!

We seriously need to parallelise this code, if it passes all the tests we subject it to.

We also need to verify the solution it gives us, and the timespan, or the simulation time for which the results are reliable and numerically stable.

The Euler Code: Solved Imbalance..but Uncertain Accuracy

The Euler method of solving differential equations is less accurate in predicting the solution than those used in the previous code. This required us to use a smaller time step of 100 seconds for time advancing.

First, I tried only the meridional flow without any diffusion or differential rotation disturbing the magnetic field. And to my surprise, I found, that not only there was no imbalance between the hemispheres, but even the flux within each hemisphere was conserved individually!

It is logically obvious, that for a symmetric initial condition, the flux distribution at any following time should be symmetric across the equator. And so it was!

Then I gave a try to the code with only differential rotation. Still no imbalance, and accurate flux conservation!

The next step was diffusion. The diffusion only code also conserved the flux individually in hemispheres, with almost no imbalance. Sometimes, there would be an imbalance, about 10 million times less than the flux in a sunspot. But then it would become zero again.

And then, I gave a run to the code, with all the 3 effects combined, and following is a screenshot showing:
1st column: Southern Hemisphere Flux
2nd column: Northern Hemisphere Flux
3rd column: Imbalance.

I also took the liberty of making a movie of the simulation, to see the flux reach the poles in 3D! The 3D plotting was one of the best things I did while in Montana.. :)

This movie was made from a simulation that runs for about 2 years in the code, and ran for about 24 hours on my laptop!

 

The video quality degraded while uploading it here. For 1080p Full HD version of the video, click here .

The Imbalance Problem

Finally, some results are here, after my workstation laptop suffered 3 reinstallations of all operating systems in 3 days!

Its good to be back after an awesome summer in Montana.

A lot of new things were tried after the last post. The major concern, however was "THE IMBALANCE" of flux which we observed.

We observed that the flux in a sunspot increases/decreases rapidly as it moves around and diffuses. The problem still hasn't been completely resolved yet.

But, to see what is wrong with the code, I tried writing a simpler version of the code following the simple Euler method of solving differential equations and simple finite difference scheme.

It was thought, that this code will be quite accurate atleast for the first few iterations, and can be used to calliberate and validate the results from the full fledged code with various numerical schemes. For simplicity, I will call the simpler code (new) the 'Euler Code' and the previous code the 'Complex Code'.

Here is an example of the problem of Imbalance:
With one sunspot each in both the hemispheres, the following screenshot shows the flux values evolving over time.
1st column: Southern Hemishpere Flux
2nd column: Northern Hemisphere Flux
3rd column: Flux imbalance between the 2 hemispheres.

Calculating the Flux Transported in different regions on the surface

Posting some results on my Birthday!

A simulation was run yesterday, with a modified profile of no. of BMRs (Unsigned Flux=2X10^21 Mx) errupting V/s time as below:

Total BMRs per cycle: 7630

The total flux between 2 latitudes was calculated by integrating the product of longitudinally averaged magnetic field and the area of a grid element at that latitude.

                              Flux=Sum[Bavg(lat)*area(lat)]

Flux between any 2 latitudes was thus calculated at regular intervals, and the difference between flux at 2 different times would give us the change in flux, meaning the flux transported. This difference can be interpreted as flux flowing in the concerned region. It can be both positive and negative.

                                                     BUTTERFLY DIAGRAM

                                                     Peak Polar Field: +/- 11 Gauss

The initial polar field in the northern hemisphere was negative, and that in the southern hemisphere was positive.

1. 0-5 :

So, the flux transported to the equator from above 5 degrees latitude would be negative in the northern hemisphere, and positive in the southern. This flux gets cancelled across the equator.

Remember that the calculation of flux transport here, shows the net flux flowing into the concerned region. So, this plot suggests that there is a flow of negative flux(leading polarity) into the northern region, and flow of positive flux(leading polarity) into the southern region. Both flow towards the equator, where these fluxes cancel each other.

The net flux in the region of 5 degrees from the equator can reach the order of 1-2 sunspots.

2. 10-20 :

The flux transport from 10-20 degree latitudes is negative in the northern hemisphere initially, and then becomes positive after half cycle. This is because the BMRs errupt at higher latitudes initially, hence the leading polarity crosses this region on its way to the equator. But, later on, BMRs start appearing at lower latitudes, which results in a net transfer of trailing polarity flux to the poles.

The net flux in this region can reach the order of  1.5X10^22 Mx.

3. 20-30 :

The trailing polarity flux from lower latitudes enters this region, signified by positive peaks in the northern hemisphere (1st cycle) and negative peaks in the southern hemisphere. 

 

The net flux in this region can reach the order of  4X10^21 Mx.

4. 30-40 :

In this region, there is just about no erruption of BMRs. So, this is like a come and go region for flux. The plot here is plotted with the unsigned difference in net flux for simplicity. Positive peaks indicate the trailing polarity flux entering the region, and negative peaks indicate, that a net flux from leading polarities has entered the region. Overall, there is a negligible storage of trailing polarity flux in this region. 

The net flux in this region can reach the order of 6X10^21 Mx. 

5. 40-50 :

There is a lot of addition and subtraction of flux in this region during a solar cycle. A positive peak in the flux transport (1st cycle) indicates the trailing polarity reaching the region, and negative peak indicates that the leading polarity flux has entered the region. Overall, there is a small net accumulation of flux from the trailing polarity in this region. 

The net flux in this region can reach the order of  8X10^21 Mx.

6. 50-60 :

This region comprises partly of the polar field. The plot shows, that there is always a positive net flux reaching this region in the first cycle, which cancels the previous field, and builds up a new one. Then, for the next cycle, there a negative flux reaching this region at all times. So, the net flux in this region does not show more than one peak. If it is increasing, it will increase till the highest value, and if it is decreasing, it will decrease till the lowest value. 

The net flux in this region can reach the order of  2.5X10^22 Mx.