Relax() stopping criterion

[jsampaio]

The relax() stopping criterion seems wrong to me. The current relax() procedure is:

1. Run until the total energy stops decreasing.
2. Run until the average torque stops decreasing.
3. Repeat step 2, decreasing MaxErr until MaxErr<=1e-9.

There are some issues that seem wrong to me:

• The torque criterion should consider the maximum torque and not the average torque. Often, only a very small subset of the spins are moving, and relax() should keep running until those spins relax. The maximum torque will measure those fastest spins, independently of the number of already relaxed spins. Moreover, the average torque is noisier, as it's a sum of lots of near-zero values.
• The torque criterion should not be that the torque stops decreasing, it should be that the torque nearly reaches zero (within the numerical noise). A slowly relaxing system may show a constant (non-zero) torque for a long time, such as a domain wall moving in a small constant field. Also, in many relaxation situations the torque may increase, before the system is relaxed. For example, a depinning domain wall at its critical field will start by deforming, slower and slower, up to a point when it depins and the torques increase -- relax() should not stop in this case.
• if the initial MaxErr is <=1e-9, the step 2 will not run.

I suggest that the maximum torque be taken instead of the average torque, and that relax() stop when maxTorque reaches a user-definable threshold. I'm not very familiar with the code, but something like changing lines 62-76 of relax.go by:

    var thresholdTorque := 1e-3; //user-settable
    maxTorque := func() float32 {
        return cuda.MaxVecNorm(solver.k1)
    }
     for maxTorque() < thresholdTorque && !pause {
            relaxSteps(N)
        }
    for MaxErr > 1e-9 && !pause {
        MaxErr /= math.Sqrt2
        for maxTorque() > thresholdTorque && !pause {
            relaxSteps(N)
        }
    }

ed: I've written the code implementing these changes in pull request #148 .

[barnex]

Although the points by @jsampaio are valid in principle, they are nearly impossible to do well in practice. The reason is that E_total and maxtorque are simply too noisy. The current implementation, although not perfect, is a least robust against noise and yields the correct results when compared to OOMMF, for a wide range of material parameters and time scales.

If desired, RunUntil(MaxTorque < myThreshold) can still be used, but this will not be robust in the sense that myThreshold will have to be chosen according to the problem. E.g.: 1e-5 might work well for permalloy, but not for a PMA material.