Same input code but different result sometimes.

ShuYuMo2003 commented 10 months ago

I run my code above serval times. But I can get different output sometimes.

I am very confuse that why I can get different output when the input code is same.

It is a bug or I set some wrong parameters?

My source code is below.

    SetGridsize(32, 8, 80)
    SetCellsize(5e-9,  5e-9, 1e-9) //Increase thickness to show the PSSW
    setPBC(5,5,1)  // cubic

    DefRegion(1,zrange(-20e-9,40e-9))
    DefRegion(2,zrange(-40e-9,-20e-9))

    //Define Regions

    Msat.SetRegion(1,196e3)   //YIG
    Aex.SetRegion(1,3.65e-12)
    alpha.SetRegion(1,6.15E-05)

    Msat.SetRegion(2,1.20E+06)    //CoFeB
    Aex.SetRegion(2,1.80E-11)
    alpha.SetRegion(2,4.00E-03)

    //Dind=0.0e-3

    // definec
    c:=1e-9
    RKKY := -1e-3 // 1mJ/m2       //  change here
    scale := (RKKY * c) / (2 * Aex.Average())
    ext_scaleExchange(1,2, scale)

    //MaxErr=1e-9

    t0 := 1/(50000000000.0)
    B_ext =vector(0.1,0.0005 * sinc(2*pi*50000000000.0*(t-t0)),0)
    minimize()

    autosave(m,1e-11)
    tableautosave(1e-11)
    tablesave()
    run(2e-09)

The output of table.txt

I got two different output. The first one is:

# t (s) mx ()   my ()   mz ()
0   0.35691237  -0.00039410926  1.0477379e-09
1.0006494964231452e-11  0.35691237  -0.0004329433   4.7086032e-05
2.0021365313454344e-11  0.3569123   -0.0003443531   -0.00013739972
3.000711253897729e-11   0.35691184  0.00013185033   -0.00036084384
4.0027933590888204e-11  0.3569123   0.00015732928   -0.00033917205
5.000563282263316e-11   0.35691318  -0.00029711242  -0.00024298255
6.000621440841988e-11   0.35691294  -0.0006298648   2.3287164e-05
7.000996931125765e-11   0.356913    -0.0004942267   0.00019759517
8.001855777044438e-11   0.3569132   -0.000104229504 0.00023100004
9.000924503796502e-11   0.35691315  8.162051e-05    0.00015210829

...

1.9000160879078468e-09  0.35691264  2.6593892e-05   -0.00015328173
1.9100027252903597e-09  0.35691223  -0.00010558424  -0.00028780213
1.9200235231559444e-09  0.35691148  -0.00011783831  -0.00033564775
1.93003744799919e-09    0.35691127  0.00010146755   -0.0003309704
1.9400280461885054e-09  0.35691148  0.00040057892   -0.00032684702
1.950017510303445e-09   0.35691163  0.0005111015    -0.0003164648
1.9600071345819544e-09  0.35691196  0.00038788197   -0.00022915583
1.9700092676709888e-09  0.3569122   0.00024653785   -5.30125e-05
1.980001081941241e-09   0.35691243  0.0002955809    0.00012328038
1.9900294383369265e-09  0.35691264  0.00046024844   0.00020312653
2e-09   0.3569123   0.0005178427    0.00020010376

I got different output after I re-run the same code.

# t (s) mx ()   my ()   mz ()
0   0.40275526  -0.57849854 -1.0178731e-12
1.0002092901318677e-11  0.4027632   -0.5785302  6.330308e-05
2.002844201567116e-11   0.40272674  -0.5784743  -0.00013651056
3.004061923046134e-11   0.40260023  -0.5781187  -0.00033430834
4.002386783718582e-11   0.40268058  -0.5780594  -0.00024187755
5.004246449862356e-11   0.4029855   -0.57836175 -0.00017168722
6.000347228811082e-11   0.40326014  -0.57864934 -2.5853999e-05
7.004030769663319e-11   0.40331417  -0.57865775 2.7038031e-05
8.00341827534107e-11    0.40317178  -0.5784423  3.0194153e-05
9.001999564526301e-11   0.40300623  -0.578272   5.2599544e-06
1.000110879746266e-10   0.40291113  -0.5782664  9.487678e-05
1.100188938439712e-10   0.4028535   -0.57834524 0.00025400924
1.2002326147317852e-10  0.40277854  -0.5784049  0.00039688707

...

1.9300117707023734e-09  0.4027431   -0.57824266 -0.00024480623
1.940026562035783e-09   0.40296608  -0.5784035  -0.00023161461
1.9500330056996213e-09  0.40318617  -0.57860154 -0.00016150001
1.9600181536443512e-09  0.40322942  -0.5785811  -5.8615384e-05
1.9700290797730503e-09  0.40310073  -0.57835704 1.7880437e-05
1.980031051139655e-09   0.40296632  -0.5781856  5.315375e-05
1.9900107231023336e-09  0.40295276  -0.57825357 9.9659526e-05
2e-09   0.40300542  -0.5784651  0.00019573218

I run my code above serval times. In most of time, the output would be the second output.

I just very confuse why I can get different output when the input code is same. It is a bug or I set some wrong parameters?

other information

I call mumax and analysis the data by python. The source code is:

import numpy as np
import scienceplots
import matplotlib.pyplot as plt
import pandas as pd
import time

np.set_printoptions(suppress=True) # Disable scientific notation for numpy
plt.style.use(['science','no-latex'])
__fig_width = 8
plt.rcParams.update({
    # 'font.family': 'sans-serif',
    # 'font.sans-serif': ['FangSong'], # , 'SimHei', 'KaiTi', 'Times New Roman', 'Tahoma', 'DejaVu Sans', 'Lucida Grande', 'Verdana'],
    # 'axes.unicode_minus': False,
    'figure.figsize': (__fig_width, __fig_width / 1.618),
    # 'legend.fontsize': 'medium',
    'axes.grid': True,
    # 'lines.linestyle': '-',
    # 'xtick.labelsize': small,
    # 'ytick.labelsize': small,
    # 'axes.labelsize': med,
    # 'figure.titlesize': large
})

def read_mumax3_table(filename):
    """Puts the mumax3 output table in a pandas dataframe"""

    from pandas import read_table

    table = read_table(filename)
    table.columns = ' '.join(table.columns).split()[1::2]

    return table

def read_mumax3_ovffiles(outputdir):
    """Load all ovffiles in outputdir into a dictionary of numpy arrays
    with the ovffilename (without extension) as key"""

    from subprocess import run, PIPE, STDOUT
    from glob import glob
    from os import path
    from numpy import load

    # convert all ovf files in the output directory to numpy files
    p = run(["mumax3-convert","-numpy",outputdir+"/*.ovf"], stdout=PIPE, stderr=STDOUT)
    if p.returncode != 0:
        print(p.stdout.decode('UTF-8'))

    # read the numpy files (the converted ovf files)
    fields = {}
    for npyfile in glob(outputdir+"/*.npy"):
        key = path.splitext(path.basename(npyfile))[0]
        fields[key] = load(npyfile)

    return fields

def run_mumax3(script, name, verbose=False):
    """ Executes a mumax3 script and convert ovf files to numpy files

    Parameters
    ----------
      script:  string containing the mumax3 input script
      name:    name of the simulation (this will be the name of the script and output dir)
      verbose: print stdout of mumax3 when it is finished
    """

    from subprocess import run, PIPE, STDOUT
    from os import path

    scriptfile = name + ".txt"
    outputdir  = name + ".out"

    # write the input script in scriptfile
    with open(scriptfile, 'w' ) as f:
        f.write(script)

    # call mumax3 to execute this script
    p = run(["mumax3","-f",scriptfile], stdout=PIPE, stderr=STDOUT)
    if verbose or p.returncode != 0:
        print(p.stdout.decode('UTF-8'))

    if path.exists(outputdir + "/table.txt"):
        table = read_mumax3_table(outputdir + "/table.txt")
    else:
        table = None

    fields = read_mumax3_ovffiles(outputdir)

    return table, fields

def detect_peaks(xx, yy, upper_x_lim, lower_y_lim):
    """
    detect peaks in the curve.
    params:
     - xx: x axis
     - yy: y axis
     - upper_x_lim: upper limit of x axis, the beyond which peaks are not detected
     - lower_y_lim: lower limit of y axix, the below which peaks are not detected
    """
    # print(list(zip(xx, yy)))
    peak_det_radius = 1
    eps_mul = 0.99
    maxidxs = []
    for idx in range(len(xx)):
        if(xx[idx] > upper_x_lim or idx - peak_det_radius < 0    \
                                or idx + peak_det_radius >= len(xx)   \
                                or yy[idx] < lower_y_lim): continue
        is_peak = True
        # is_peak = yy[idx] > yy[idx - 1] and yy[idx] > yy[idx + 1]
        for t in range(peak_det_radius):
            if(yy[idx - t - 1] > yy[idx - t] * eps_mul):
                is_peak = False
        for t in range(peak_det_radius):
            if(yy[idx + t + 1] > yy[idx + t] * eps_mul):
                is_peak = False
        if(is_peak):
            maxidxs.append(idx)
    pts = [(xx[i], yy[i]) for i in maxidxs]
    pts.sort(key=lambda x : x[1], reverse=True)
    return pts

def merge_ticks(original_xticks, extra_xticks, raduis=5):
    ticks = extra_xticks.copy()
    for tick in original_xticks:
        can_be_added = True
        for orj in ticks:
            if(abs(tick - orj) < raduis):
                can_be_added = False
                break
        if(can_be_added):
            ticks.append(tick)
    return ticks

def run_specific_mumax_script(fmax, T, dt, uniform_name):
    # Note that this is a format string, this means that the statements inside the
    # curly brackets get evaluated by python. In this way, we insert the values of
    # the variables above in the script.
    script=f"""

    // Set solver type. 1:Euler, 2:Heun, 3:Bogaki-Shampine, 4: Runge-Kutta (RK45), 5: Dormand-Prince, 6: Fehlberg, -1: Backward Euler
    // SetSolver(1) 

    SetGridsize(32, 8, 80)
    SetCellsize(5e-9,  5e-9, 1e-9) //Increase thickness to show the PSSW
    setPBC(5,5,1)  // cubic

    DefRegion(1,zrange(-20e-9,40e-9))
    DefRegion(2,zrange(-40e-9,-20e-9))

    //Define Regions

    Msat.SetRegion(1,196e3)   //YIG
    Aex.SetRegion(1,3.65e-12)
    alpha.SetRegion(1,6.15E-05)

    Msat.SetRegion(2,1.20E+06)    //CoFeB
    Aex.SetRegion(2,1.80E-11)
    alpha.SetRegion(2,4.00E-03)

    //Dind=0.0e-3

    // definec
    c:=1e-9
    RKKY := -1e-3 // 1mJ/m2       //  change here
    scale := (RKKY * c) / (2 * Aex.Average())
    ext_scaleExchange(1,2, scale)

    //MaxErr=1e-9

    t0 := 1/({fmax})
    B_ext =vector(0.1,0.0005 * sinc(2*pi*{fmax}*(t-t0)),0)
    minimize()

    autosave(m,{dt})
    tableautosave({dt})
    tablesave()
    run({T})
    """

    table, fields = run_mumax3(script,uniform_name)
    return table, fields

def main():
    uniform_name = 'CACHE_DATA\\' + 'TEST_' + time.strftime("%Y-%m-%d_%H-%M-%S", time.localtime())

    # NUMERICAL PARAMETERS RELEVANT FOR THE SPECTRUM ANALYSIS
    fmax = 50e9        # maximum frequency (in Hz) of the sinc pulse
    T    = 2e-9        # simulation time (longer -> better frequency resolution)
    dt   = 1/(2*fmax)  # the sample time (Nyquist theorem taken into account)

    print(uniform_name, 'running script.')
    table, fields = run_specific_mumax_script(fmax, T, dt, uniform_name)

    print(uniform_name, 'running fft.')
    # FAST FOURIER TRANSFORM
    dm     = table["mz"] - table["mz"][0]   # average magnetization deviaton
    spectr = np.abs(np.fft.fft(dm))         # the absolute value of the FFT of dm
    freq   = np.linspace(0, 1/dt, len(dm))  # the frequencies for this FFT

    print(uniform_name, 'plotting.')
    # PLOT THE SPECTRUM
    xx = freq / 1e9
    yy = spectr
    upper_x_lim = fmax/1e9

    fig, ax = plt.subplots()

    ax.plot(xx, yy)

    peaks = detect_peaks(xx, yy, upper_x_lim, 0)
    ylim_before = ax.get_ylim()
    for peak in peaks:
        ax.plot([peak[0], peak[0]], [-5, peak[1]], color="k",linestyle=':')
    ax.set_xlim(0,upper_x_lim)
    ax.set_ylim(ylim_before)

    original_xticks = ax.get_xticks()
    extra_xticks = list(map(lambda x : x[0], peaks))
    ax.set_xticks(merge_ticks(original_xticks, extra_xticks, 4))

    ax.set_ylabel("Spectrum (a.u.)")
    ax.set_xlabel("Frequency (GHz)")
    plt.xticks(rotation=46) # rotate xticks
    print('peaks:')
    pd.DataFrame(sorted(peaks), columns=["Frequency (GHz)", "Spectrum (a.u.)"], index=range(1, len(peaks) + 1)).to_csv(uniform_name + '-peaks.csv')
    plt.savefig(uniform_name + '.png')
    pd.DataFrame(list(zip(xx, yy))).to_csv(uniform_name + '-xy.csv')

    print(uniform_name, 'process done.')

while(True):
    main()

My environment is: Processor Intel(R) Xeon(R) W-2225 CPU @ 4.10GHz 4.10 GHz Installed RAM 64.0 GB (63.6 GB usable) System type 64-bit operating system, x64-based processor GPU Quadro P2200(5119MB), CUDA Driver 12.2, cc=6.1, using cc=61 PTX

MathieuMoalic commented 10 months ago

Hi,

Mumax is not deterministic in that way, I'm almost sure it's because of CPU/GPU floating point precision mismatch. This leads to tiny variations of the starting magnetizations between runs which can lead to huge changes if you have an unstable system. A workaround is to save the starting magnetization and load it, the loaded magnetization will always lead to the same outcome of the similation. Let's hope this won't be overlooked in the upcoming mumax4 release.

ShuYuMo2003 commented 10 months ago

Hiii, Thanks for your explaination !

I will learn and try to save the initial magnetization and reload it.

But I still confuse that the floating point precision mismatch may always lead the same mistake for the same input. This may not explain the different output in this system. I guass this behaviour may be caused by some randomness in the mumax?

jplauzie commented 10 months ago

Hi,

For questions about scripts, it's best to post them in the mumax group here (it is named after mumax2, but is still active for mumax3).

That said: mumax is indeed deterministic up to single precision. So using the same input will always give the same output, unless your system is sensitive. There are a few psuedo-RNGs, but you can control the seeds using variables like randSeed, RandomMagSeed, and ThermSeed. In your case, these aren't affecting your simulation, just mentioning for completeness.

Very occasionally, you will have problems that are sensitive to the precision. It's usually pretty rare, however. I think it usually happens when there is a symmetry that isn't broken. In those cases, it's like a marble on top of a hill- it's an unstable equilibrium, and the noise from single precision is enough to get it rolling one way or the other. You should probably save and compare the total energies to see which is the true ground state.

Can you give some more details, like what version of mumax are you using, and how often you see the different results? And what the magnetization states look like? Guessing from just the average magnetizations, I would suspect the (0.35691237, -0.00039410926,1.0477379e-09) state is probably two layers are anti-parallel. So probably minimize can't decide if they should be anti-parallel or parallel, or something.

Looking at your output, it seems like there are different initial states in the table even before run(). Probably minimize() is having trouble deciding which state is better (you're also not setting an initial magnetization, so you're starting from a high energy state. You would probably be better of starting with relax() first, or initializing something close to the expected state. There is a warning in the Examples page about this. You might still run into issues, relax still has sensitivity to numerical rounding but it might help). Also/alternatively, you can probably nudge it towards the proper ground state by applying a small static field component in the y-direction (keeping in mind that minimize/relax disable excitations, so the sinc you provide is disabled when relax/minimize are running). If either of those two things doesn't solve it, there are some other tricks like adding some noise you can do to try to nudge it towards the correct state.

-Separately, you might want to fully simulate your system, instead of using PBC. But that I am less sure about.

Also, you can look to issues like this one, which discuss similar things: https://github.com/mumax/3/issues/150 https://github.com/mumax/3/pull/233

Numerical rounding is a known issue, and it's further compounded since GPU tasks aren't always deterministic, either.

As for mumax using both double/single precision: It's true that some of the CPU calculations are done in double precision (and then truncated), and the GPU are done in single. But I think this was an intentional design decision, and hard to avoid. Single precision is necessary to optimize mumax for performance on consumer-grade GPUs. Unfortunately decent double-precision performance is locked behind workstation tier cards (with a few discontinued exceptions).

Best, Josh L.

ShuYuMo2003 commented 10 months ago

@jplauzie Very thanks for your detailed explaination and suggestions.

I set the a proper initial magnetization. And mumax gave the reasonable output.

Again, thank you for your help !!

mumax / 3