0

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:

  1. 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.
  2. 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!

  • [Way to go from recursion to iteration](https://stackoverflow.com/q/159590/327083) – J... Jan 12 '22 at 17:35
  • [Can every recursion be converted into iteration?](https://stackoverflow.com/q/931762/327083) – J... Jan 12 '22 at 17:36
  • Recursion is a bad idea, almost by definition, when it results in a solution that explodes like this. Avoiding it will be your solution. A rewrite probably isn't what you were hoping for, but recursive algorithms really should come with disclaimers and warnings - they can be effective in certain cases, but they put a pretty hard limit on the complexity of problem you can solve, so it's easy to paint yourself into a corner with a bad early decision to go that way. – J... Jan 12 '22 at 17:38
  • Oh this is actually perfectly helpful. It is quite easy to rewrite. Thanks a lot! – Jan Mulder Jan 12 '22 at 18:00

0 Answers0