Page MenuHomec4science
Files F91676826

simulation.cc
No OneTemporary
Actions

File Metadata

Created
Wed, Nov 13, 09:06

simulation.cc
View Options

/* -------------------------------------------------------------------------- */
#include "simulation.hh"
/* -------------------------------------------------------------------------- */
#include <cmath>
#include <iostream>
#include <omp.h>
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
Simulation::Simulation(int m, int n)
: m_global_m(m), m_global_n(n), m_epsilon(1e-7), m_h_m(1. / m),
m_h_n(1. / n), m_grids(m, n), m_f(m, n),
m_dumper(new DumperASCII(m_grids.old())) {
m_n_threads = omp_get_max_threads();
}
/* -------------------------------------------------------------------------- */
void Simulation::set_initial_conditions() {
#pragma omp for schedule(static, m_global_m/m_n_threads)
for (int i = 0; i < m_global_m; i++) {
for (int j = 0; j < m_global_n; j++) {
m_f(i, j) = -2. * 100. * M_PI * M_PI * std::sin(10. * M_PI * i * m_h_m) *
std::sin(10. * M_PI * j * m_h_n);
}
}
}
/* -------------------------------------------------------------------------- */
std::tuple<float, int> Simulation::compute() {
int s = 0;
do {
#pragma omp barrier
#pragma omp single
l2 = 0.;
compute_step();
#pragma omp single nowait
m_grids.swap();
++s;
} while (l2 > m_epsilon);
#pragma omp single nowait
m_dumper->dump(s);
return std::make_tuple(std::sqrt(l2), s);
}
/* -------------------------------------------------------------------------- */
void Simulation::set_epsilon(float epsilon) { m_epsilon = epsilon * epsilon; }
/* -------------------------------------------------------------------------- */
float Simulation::epsilon() const { return std::sqrt(m_epsilon); }
/* -------------------------------------------------------------------------- */
void Simulation::compute_step() {
Grid & u = m_grids.current();
Grid & uo = m_grids.old();
#pragma omp for reduction(+:l2) schedule(static, m_global_m/m_n_threads)
for (int i = 1; i < m_global_m - 1; i++) {
for (int j = 1; j < m_global_m - 1; j++) {
// computation of the new step
u(i, j) = 0.25 * (uo(i - 1, j) + uo(i + 1, j) + uo(i, j - 1) +
uo(i, j + 1) - m_f(i, j) * m_h_m * m_h_n);
// L2 norm
l2 += (uo(i, j) - u(i, j)) * (uo(i, j) - u(i, j));
}
}
}

Event Timeline

c4science · Help