diff options
author | Jaron Kent-Dobias <jaron@kent-dobias.com> | 2018-10-15 22:57:17 -0400 |
---|---|---|
committer | Jaron Kent-Dobias <jaron@kent-dobias.com> | 2018-10-15 22:57:17 -0400 |
commit | 1343a3fe6bd17a2487f12a0d61be8dc83cd722a0 (patch) | |
tree | fa937f0f3ba0f4977036c862c846a2ee461540ca /examples/simple_ising.cpp | |
parent | 6e8b19e1f1a244ef09e1b63d7593250d6ce01692 (diff) | |
download | c++-1343a3fe6bd17a2487f12a0d61be8dc83cd722a0.tar.gz c++-1343a3fe6bd17a2487f12a0d61be8dc83cd722a0.tar.bz2 c++-1343a3fe6bd17a2487f12a0d61be8dc83cd722a0.zip |
many changes, including reworking the measurements system
Diffstat (limited to 'examples/simple_ising.cpp')
-rw-r--r-- | examples/simple_ising.cpp | 192 |
1 files changed, 192 insertions, 0 deletions
diff --git a/examples/simple_ising.cpp b/examples/simple_ising.cpp new file mode 100644 index 0000000..24e4ae5 --- /dev/null +++ b/examples/simple_ising.cpp @@ -0,0 +1,192 @@ + +#include <getopt.h> +#include <iostream> + +#include "include/randutils/randutils.hpp" + +#include <wolff.hpp> + +// define your R_t and X_t. here both are the same, ising_t +class ising_t { + public: + bool x; + + // both X_t and R_t need a default constructor (and destructor, if relevant) + ising_t() : x(false) {} + ising_t(bool x) : x(x) {} + + // R_t needs the member functions + // X_t act(const X_t& s) const {} + // R_t act(const R_t& s) const {} + // to define action on both spins and other transformations + ising_t act(const ising_t& s) const { + if (x) { + return ising_t(!(s.x)); + } else { + return ising_t(s.x); + } + } + + // R_t needs the member functions + // X_t act_inverse(const X_t& s) const {} + // R_t act_inverse(const R_t& s) const {} + // to define action of its inverse on both spins and other transformations + ising_t act_inverse(const ising_t& s) const { + return this->act(s); + } +}; + +// define how measurements should be taken by importing the interface wolff_measurement<R_t, X_t> +class ising_measurements : public wolff_measurement<ising_t, ising_t> { + private: + count_t n; + + double E; + int M; + v_t S; + + double totalE; + double totalM; + double totalS; + + public: + ising_measurements(D_t D, L_t L, double H) { + n = 0; + M = (int)pow(L, D); + E = -D * pow(L, D) - H * pow(L, D); + + totalE = 0; + totalM = 0; + totalS = 0; + } + + void pre_cluster(const state_t<ising_t, ising_t>& s, count_t step, count_t N, v_t v0, const ising_t& R) { + S = 0; + } + + void plain_bond_added(v_t v, const ising_t& s_old, const ising_t& s_new, v_t vn, const ising_t& sn, double dE) { + E += dE; + } + + void ghost_bond_added(v_t v, const ising_t& s_old, const ising_t& s_new, double dE) { + E += dE; + + if (s_old.x) { + M++; + } else { + M--; + } + + if (s_new.x) { + M--; + } else { + M++; + } + } + + void plain_site_transformed(v_t v, const ising_t& s_old, const ising_t& s_new) { + S++; + } + + void ghost_site_transformed(const ising_t& R_old, const ising_t& R_new) { + } + + void post_cluster(const state_t<ising_t, ising_t>& s, count_t step, count_t N) { + totalE += E; + totalM += M; + totalS += S; + n++; + } + + double avgE() { + return totalE / n; + } + + double avgM() { + return totalM / n; + } + + double avgS() { + return totalS / n; + } +}; + +int main(int argc, char *argv[]) { + + // set defaults + count_t N = (count_t)1e4; + D_t D = 2; + L_t L = 128; + double T = 2.26918531421; + double H = 0.0; + + int opt; + + // take command line arguments + while ((opt = getopt(argc, argv, "N:D:L:T:H:")) != -1) { + switch (opt) { + case 'N': // number of steps + N = (count_t)atof(optarg); + break; + case 'D': // dimension + D = atoi(optarg); + break; + case 'L': // linear size + L = atoi(optarg); + break; + case 'T': // temperature + T = atof(optarg); + break; + case 'H': // external field + H = atof(optarg); + break; + default: + exit(EXIT_FAILURE); + } + } + + // define the spin-spin coupling + std::function <double(const ising_t&, const ising_t&)> Z = [] (const ising_t& s1, const ising_t& s2) -> double { + if (s1.x == s2.x) { + return 1.0; + } else { + return -1.0; + } + }; + + // define the spin-field coupling + std::function <double(const ising_t&)> B = [=] (const ising_t& s) -> double { + if (s.x) { + return -H; + } else { + return H; + } + }; + + // initialize the system + state_t<ising_t, ising_t> s(D, L, T, Z, B); + + // initialize the random number generator + randutils::auto_seed_128 seeds; + std::mt19937 rng{seeds}; + + // define function that generates self-inverse rotations + std::function <ising_t(std::mt19937&, const ising_t&)> gen_R = [] (std::mt19937&, const ising_t& s) -> ising_t { + return ising_t(true); + }; + + // initailze the measurement object + ising_measurements m(D, L, H); + + // run wolff N times + wolff<ising_t, ising_t>(N, s, gen_R, m, rng); + + // print the result of our measurements + std::cout << "Wolff complete!\nThe average energy per site was " << m.avgE() / s.nv + << ".\nThe average magnetization per site was " << m.avgM() / s.nv + << ".\nThe average cluster size per site was " << m.avgS() / s.nv << ".\n"; + + // exit + return 0; +} + |