I am applying the so-called Wolff algorithm to simulate a very large square lattice Ising model at critical temperature in Visual Studio C#. I start out with all lattice sites randomly filled with either -1 or +1, and then the Wolff algorithm is applied for a certain number of times. Basically, per application, it randomly activates one lattice site and it checks if neighbouring lattice sites are of the same spin (the same value as the initially selected site) and with a certain chance it gets activated as well. The process is repeated for every new activation, so the activated lattice sites form a connected subset of the full lattice. All activated spins flip (*= -1). This is done a large number of times.
This is the class I use:
class WolffAlgorithm
{
public sbyte[,] Data;
public int DimX;
public int DimY;
public double ExpMin2BetaJ;
private sbyte[,] transl = { { -1, 0 }, { 1, 0 }, { 0, -1 }, { 0, 1 } };
public WolffAlgorithm(sbyte[,] data, double beta, double j)
{
Data = data;
DimX = data.GetLength(0);
DimY = data.GetLength(1);
ExpMin2BetaJ = Math.Exp(- 2 * beta * j);
}
public void Run(int n, int print_n)
{
Random random = new Random();
DateTime t0 = DateTime.Now;
int act = 0;
for (int i = 0; i < n; i++)
{
int x = random.Next(0, DimX);
int y = random.Next(0, DimY);
int[] pos = { x, y };
List<int[]> activations = new List<int[]>();
activations.Add(pos);
sbyte up_down = Data[x, y];
Data[x, y] *= -1;
CheckActivation(random, up_down, activations, x, y);
act += activations.Count;
if ((i + 1) % print_n == 0)
{
DateTime t1 = DateTime.Now;
Console.WriteLine("n: " + i + ", act: " + act + ", time: " + (t1 - t0).TotalSeconds + " sec");
t0 = t1;
act = 0;
}
}
}
private void CheckActivation(Random random, sbyte up_down, List<int[]> activations, int x, int y)
{
for (int i = 0; i < transl.GetLength(0); i++)
{
int tr_x = (x + transl[i, 0]) % DimX;
if (tr_x < 0) tr_x += DimX;
int tr_y = (y + transl[i, 1]) % DimY;
if (tr_y < 0) tr_y += DimY;
if (Data[tr_x, tr_y] != up_down)
{
continue;
}
if (random.NextDouble() < 1 - ExpMin2BetaJ)
{
int[] new_activation = { tr_x, tr_y };
activations.Add(new_activation);
Data[tr_x, tr_y] *= -1;
CheckActivation(random, up_down, activations, tr_x, tr_y);
}
}
}
}
When I approach the equilibrium configuration, the activated groups become so large, that the recursive function exceeds the maximum number of stack frames (apparently 5000), triggering the StackOverflowException (fitting for my first post on StackOverflow.com, I guess haha).
What I've tried:
- I have tried adjusting the maximum number of stack frames like what is proposed here: https://www.poppastring.com/blog/adjusting-stackframes-in-visual-studio This did not work for me.
- I have also tried to reduce the number of stack frames needed by just copying the same function in the function (however ugly this looks):
private void CheckActivation2(Random random, sbyte up_down, List<int[]> activations, int x, int y)
{
for (int i = 0; i < transl.GetLength(0); i++)
{
int tr_x = (x + transl[i, 0]) % DimX;
if (tr_x < 0) tr_x += DimX;
int tr_y = (y + transl[i, 1]) % DimY;
if (tr_y < 0) tr_y += DimY;
if (Data[tr_x, tr_y] != up_down)
{
continue;
}
if (random.NextDouble() < 1 - ExpMin2BetaJ)
{
int[] new_activation = { tr_x, tr_y };
activations.Add(new_activation);
Data[tr_x, tr_y] *= -1;
for (int i2 = 0; i2 < transl.GetLength(0); i2++)
{
int tr_x2 = (tr_x + transl[i2, 0]) % DimX;
if (tr_x2 < 0) tr_x2 += DimX;
int tr_y2 = (tr_y + transl[i2, 1]) % DimY;
if (tr_y2 < 0) tr_y2 += DimY;
if (Data[tr_x2, tr_y2] != up_down)
{
continue;
}
if (random.NextDouble() < 1 - ExpMin2BetaJ)
{
int[] new_activation2 = { tr_x2, tr_y2 };
activations.Add(new_activation2);
Data[tr_x2, tr_y2] *= -1;
CheckActivation2(random, up_down, activations, tr_x2, tr_y2);
}
}
}
}
}
This works, but it is not enough and I don't know how to scale this elegantly by an arbitrary number of times, without calling the function recursively.
I hope anyone can help me with this!