Projekt 2: Monte-Carlo Simulation

21.06.2019 / B. Leder

Wissenschaftliches Rechnen III / CP III

Magnetisierung

Bestimmen Sie \(|M_0|^2=<|M|^2>\) bei \(h=0\) und \(\lambda=1.0\) für \(d=3\), \(L=4, 8, 16\) und jeweils mindestens 20 Werte für \(\kappa\). Wählen Sie die \(\kappa\)-Werte um den Meanfield-Wert \(\approx 0.2333\) für das kritische \(\kappa_c\). Interpretieren Sie Ihre Daten (Abhängigkeit von \(L\) und \(\kappa\)) und geben Sie einen Schätzwert für \(\kappa_c\) an.

Anforderungen / Hinweise:

#include <cuComplex.h>
...
cuDoubleComplex *phi;
phi=(cuDoubleComplex*)malloc(nvol*sizeof(cuDoubleComplex));

Zur Verwendung von cuDoubleComplex sehen Sie sich /usr/local/cuda/include/cuComplex.h an. Der fortlaufende Index sei gegeben als \(i=x_1+N_1\,x_2+\dots+(N_1\,N_2\dots N_d)x_d\). Zur Bestimmung des Index des nächsten Nachbarn mit periodischen Randbedingungen können sie geom_pbc.cu verwenden, siehe bereitgestelltes Archiv [1].

#include "stat5.h"
...
// Initialisierung des Statistikmoduls, eine Observable
clear5(1,500);
// Sweeps
for (i=1; i<=nsweep; i++)
{
  metro_sweep(delta);
  m=magnet();
  // Messwert für Auswertung speichern
  accum5(1,cuCabs(m)*cuCabs(m));
}
// Mittelwert, Fehler, Auto-Korrelationszeit
printf(" |M|^2:    %f, %f, %f\n",aver5(1),sigma5(1),tauint5(1));

Monte-Carlo-Simulation: GPU

Planen Sie die Implementierung auf der GPU. Dabei soll ein Sweep und die Messung vollständig auf der GPU laufen. Das Ergebnis einer Messung soll sofort auf die CPU kopiert werden. Geben Sie eine Liste der notwendigen Kernel an und beschreiben Sie jeweils kurz die Aufgaben und die geplante Threadverteilung (Datenverteilung auf Threads). Zeichnen Sie ein Programmablaufplan.

Beachten Sie hierbei besonders, dass der Sweep die Gleichgewichtsbedingung erfüllt. Im Allgemeinen hängen die Updates der einzelnen Spins voneinander ab (z.B. wenn einer der Nachbarn bereits ein Update erfahren hat). Dann muss die Reihenfolge in jedem Sweep gleich sein. Das ist ungünstig für die Implementierung auf der GPU (parallele Ausführung). Ein Ausweg ist die Punkte des Gitters in Untermengen aufzuteilen, so dass die Updates der Spins in einer Untermenge unabhängig voneinander sind. Diese können dann in beliebiger Reihenfolge durchlaufen werden und nur die Abfolge der Untermengen muss für jeden Sweep gleich sein.

Implementieren Sie die Aktualiserung aller Spins (ein sweep) auf der GPU. Die Messung der Observablen nach jedem sweep (oder einer festen Anzahl) kann auf der CPU erfolgen. So ist für große Gitter ein speedup von \(>10\) möglich.

Anforderungen / Hinweise

int *d_nn;
CHECK(cudaMalloc((void**)&d_nn,nvol*(2*ndim+1)*sizeof(int)));
CHECK(cudaMemcpy(d_nn, nn[0], nvol*(2*ndim+1)*sizeof(int), cudaMemcpyHostToDevice));
// ausserhalb von Funktion/Kernel
__device__ double devLambda, devKappa;
...
// in einer Host-Funktion
CHECK(cudaMemcpyToSymbol(devLambda, &lambda, sizeof(double)));
CHECK(cudaMemcpyToSymbol(devKappa, &kappa, sizeof(double)));

devLambda und devKappa können dann in jedem Kernel verwendet werden.

Monte-Carlo-Simulation: Phasen in \(d=2,3,4\)

Die unterschiedlichen Phasen lassen sich durch das Verhalten der Magnetisierung in Abhängigkeit von der Systemgröße \(L^d\) (finite size scaling) kennzeichnen:

\[<|M|^2>\overset{L\to\infty}{=} \begin{cases} \chi\, L^{-d} & \text{ungeordnet/paramagnetisch/symmetrisch}\\ \text{const}\, L^{-\eta} & \text{Kosterlitz–Thouless} \\ |M_0|^2 & \text{ferromagnetisch/Goldstone} \end{cases}\]

mit der Suszeptibilität \(\chi\) und \(\eta=O(\frac{1}{\kappa}\)) (im XY-Model \(\eta=\frac{1}{4\pi\kappa}\)).

Die Kosterlitz–Thouless Phase tritt in \(d=2\) auf (siehe auch Physik-Nobelpreis 2016).

Zur Untersuchung der Phasen definieren wir das Verhältnis der Magnetisierung bei \(2L\) und \(L\)

\[R(L)=\frac{<|M|^2>_{2L}}{<|M|^2>_{L}}\]

Aufgaben:

  1. Leiten Sie für die unterschiedlichen Phasen das Verhalten von \(R(L)\) für \(L\to\infty\) ab.

  2. In MC Simulationen bestimmen Sie \(<|M|^2>\) als Mittelwert mit einem endlichen statistischen Fehler. Warum können Sie einfache Fehlerfortpflanzung für den Quotienten \(R(L)\) benutzen (im Gegensatz z.B. zum Binder-Parameter)?

  3. Bestimmen Sie (mit Hilfe von \(R(L)\)) \(\kappa_c\) für \(\lambda=1.0\) und \(d=2,3,4\), d.h. das kritische \(\kappa\), das die paramagnetische Phase von der ferromagnetischen Phase trennt. Bestimmen Sie dafür den Schnittpunkt von \(R(L_1)\) und \(R(L_2)\).

  4. Verifizieren Sie (im Rahmen der statistichen Fehler) die Abhängikeit von \(R(\infty)\)

Hinweis: Sie sollten \(R(L)\) jeweils für mindestens 3 Werte von \(L\) bestimmen.

[1] Vorgegebene CPU-Version http://people.physik.hu-berlin.de/~leder/cp3/uebung8.zip