From 53f7c5b3fe72249d3fca2f58d1d2637a7f6f3464 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=D0=9F=D0=B8=D1=81=D0=B0=D1=80=D0=B5=D0=B2=20=D0=92=D0=B0?=
=?UTF-8?q?=D1=81=D0=B8=D0=BB=D0=B8=D0=B9=20=D0=92=D1=8F=D1=87=D0=B5=D1=81?=
=?UTF-8?q?=D0=BB=D0=B0=D0=B2=D0=BE=D0=B2=D0=B8=D1=87?=
<vpisarev@miem.hse.ru>
Date: Wed, 6 Oct 2021 12:55:03 +0000
Subject: [PATCH] Added 1D heat equation solvers
---
heat_equation.hpp | 60 ++++++++++
heat_equation_cg.cpp | 208 +++++++++++++++++++++++++++++++++
heat_equation_explicit.cpp | 156 +++++++++++++++++++++++++
heat_equation_explicit_omp.cpp | 157 +++++++++++++++++++++++++
heat_equation_implicit.cpp | 188 +++++++++++++++++++++++++++++
5 files changed, 769 insertions(+)
create mode 100644 heat_equation.hpp
create mode 100644 heat_equation_cg.cpp
create mode 100644 heat_equation_explicit.cpp
create mode 100644 heat_equation_explicit_omp.cpp
create mode 100644 heat_equation_implicit.cpp
diff --git a/heat_equation.hpp b/heat_equation.hpp
new file mode 100644
index 0000000..350efc9
--- /dev/null
+++ b/heat_equation.hpp
@@ -0,0 +1,60 @@
+#ifndef HEAT_EQN_H
+#define HEAT_EQN_H
+
+#include "../vectors_and_matrices/array_types.hpp"
+#include <iostream>
+#include <fstream>
+#include <cstdint>
+
+using intptr_t = std::intptr_t;
+
+// Общее граничное условие вида
+// (a u + b du/dx)_bound = z
+struct boundary_cond
+{
+ double coeff_u;
+ double coeff_du;
+ double rhp_val;
+};
+
+struct heat_eq_problem
+{
+ vec<double> us;
+ vec<double> us_next;
+ double step_x;
+ double conductivity;
+ boundary_cond left_bc;
+ boundary_cond right_bc;
+};
+
+boundary_cond dirichlet_bc(double val)
+{
+ return boundary_cond {1.0, 0.0, val};
+}
+
+boundary_cond neumann_bc(double val)
+{
+ return boundary_cond {0.0, 1.0, val};
+}
+
+void print_temp_profile(std::ostream& log_file, heat_eq_problem prob)
+{
+ double h = prob.step_x;
+ intptr_t nx = prob.us.length() - 2;
+ log_file << "0.0 " << 0.5 * (prob.us(0) + prob.us(1)) << "\n";
+ for (intptr_t i = 1; i <= nx; i++)
+ {
+ log_file << (i - 0.5) * h << " " << prob.us(i) << "\n";
+ }
+ log_file << nx * h << " " << 0.5 * (prob.us(nx) + prob.us(nx+1)) << "\n";
+ log_file << "\n\n" << std::endl;
+}
+
+void print_log(std::ostream& log_file, intptr_t step, double step_time,
+ heat_eq_problem prob)
+{
+ log_file << "Timestep " << step << ". t = " << step_time << "\n";
+ print_temp_profile(log_file, prob);
+}
+
+#endif
diff --git a/heat_equation_cg.cpp b/heat_equation_cg.cpp
new file mode 100644
index 0000000..b4356f0
--- /dev/null
+++ b/heat_equation_cg.cpp
@@ -0,0 +1,208 @@
+#include "../vectors_and_matrices/array_types.hpp"
+#include "../vectors_and_matrices/linalg.hpp"
+#include "heat_equation.hpp"
+#include <iostream>
+#include <fstream>
+#include <cstddef>
+#include <string>
+#include <cmath>
+#include <omp.h>
+
+using size_t = std::size_t;
+using intptr_t = std::intptr_t;
+
+void set_constant_u(vec<double> us, double u0)
+{
+ intptr_t n = us.length()-2;
+ for (intptr_t i=1; i<=n; i++)
+ {
+ us(i) = u0;
+ }
+}
+
+void compute_residual(vec<double> residual, vec<double> x,
+ vec<double> lower, vec<double> diag, vec<double> upper,
+ vec<double> rhp
+ )
+{
+ intptr_t n = rhp.length() - 2;
+ for (intptr_t i=1; i<=n; i++)
+ {
+ residual(i) = lower(i) * x(i-1) + diag(i) * x(i) + upper(i) * x(i+1) - rhp(i);
+ }
+}
+
+double pAp(vec<double> dir,
+ vec<double> lower, vec<double> diag, vec<double> upper
+ )
+{
+ // Вычисляет dir' * A * dir,
+ // где A - трёхдиагональная матрица, задаваемая l, d и u
+ double prod = 0;
+ intptr_t n = dir.length() - 2;
+ for (intptr_t i=1; i<=n; i++)
+ {
+ prod += dir(i) * (lower(i) * dir(i-1) + diag(i) * dir(i) + upper(i) * dir(i+1));
+ }
+ return prod;
+}
+
+void solve_tridiag(vec<double> x,
+ vec<double> lower, vec<double> diag, vec<double> upper,
+ vec<double> rhp
+ )
+{
+ // входные параметры:
+ // lower, diag, upper — массивы коэффициентов l, d, u
+ // rhp — массив с правой частью (right-hand part)
+ // x — массив, в который будет записываться решение
+ // courant - a * tau / h^2
+
+ // РїСЂСЏРјРѕР№ С…РѕРґ
+ double factor;
+ intptr_t n = rhp.length()-2;
+ vec<double> residual(n+2);
+ vec<double> dir(n+2);
+
+ diag(1) -= upper(0) * lower(1) / diag(0);
+ rhp(1) -= rhp(0) * lower(1) / diag(0);
+ diag(n) -= lower(n+1) * upper(n) / diag(n+1);
+ rhp(n) -= rhp(n+1) * upper(n) / diag(n+1);
+
+ lower(1) = 0;
+ upper(n) = 0;
+
+ compute_residual(residual, x, lower, diag, upper, rhp);
+
+ residual(0) = 0;
+ residual(n+1) = 0;
+ for (intptr_t i = 0; i <= n+1; i++)
+ {
+ dir(i) = -residual(i);
+ }
+
+ double rkrk = dot(residual, residual);
+ while (sqrt(rkrk / n) > 1e-6)
+ {
+ double alpha = rkrk / pAp(dir, lower, diag, upper);
+
+ // x <- x + alpha * dir
+ // r <- r + alpha * A * dir
+ for (intptr_t i = 1; i <= n; i++)
+ {
+ x(i) += alpha * dir(i);
+ residual(i) += alpha * (lower(i) * dir(i-1) + diag(i) * dir(i) + upper(i) * dir(i+1));
+ }
+
+ double d = dot(residual, residual);
+ double beta = d / rkrk;
+ rkrk = d;
+
+ // dir <- -r + beta * dir
+ for (intptr_t i = 1; i <= n; i++)
+ {
+ dir(i) = -residual(i) + beta * dir(i);
+ }
+ }
+}
+
+void dirichlet_implicit_step(heat_eq_problem& prob, double tau,
+ double u0, double uL,
+ vec<double> lower, vec<double> diag, vec<double> upper)
+{
+ vec<double> us = prob.us;
+ intptr_t n = us.length() - 2;
+ double h = prob.step_x;
+ double courant = tau * prob.coeff_a / (h * h);
+ for (intptr_t i = 0; i <= n; i++)
+ {
+ diag(i) = 1.0 + 2 * courant;
+ lower(i) = -courant;
+ upper(i) = -courant;
+ }
+ upper(0) = 1.0;
+ lower(n+1) = 1.0;
+ diag(0) = 1.0;
+ diag(n+1) = 1.0;
+ us(0) = 2 * u0;
+ us(n+1) = 2 * uL;
+ solve_tridiag(prob.us_next, lower, diag, upper, us);
+
+ vec<double> tmp = us;
+ prob.us = prob.us_next;
+ prob.us_next = tmp;
+ apply_dirichlet_bc(prob, u0, uL);
+}
+
+void dirichlet_implicit_solver(heat_eq_problem& prob, double u0, double uL,
+ double L, double tau, size_t ntimesteps,
+ size_t log_freq, std::ofstream& log_file)
+{
+ // входные параметры:
+ // prob - начальная конфигурация
+ // u0 - температура при x=0
+ // uL - температура при x=L
+ // L - длина отрезка, на котором решается уравнение
+ // tau - шаг по времени
+ // ntimesteps - сколько шагов по времени сделать
+ // log_freq - раз в столько шагов выводить информацию
+ // log_file - файловый дескриптор для вывода информации
+ size_t nx = prob.us.length() - 2;
+ double h = L / nx;
+ prob.step_x = h;
+ double courant = prob.coeff_a * tau / (h * h);
+ if (courant > 0.5) {
+ std::cerr << "Warning: Courant number > 1\n";
+ }
+ print_log(log_file, 0, 0.0, prob);
+ vec<double> lower(nx+2);
+ vec<double> diag(nx+2);
+ vec<double> upper(nx+2);
+ set_constant_u(lower, -courant);
+ set_constant_u(upper, -courant);
+ upper(0) = 1.0;
+ lower(nx+1) = 1.0;
+ // эволюция по времени
+ for (intptr_t step = 1; step <= ntimesteps; step++) {
+ dirichlet_implicit_step(prob, tau, u0, uL, lower, diag, upper);
+ if (step % log_freq == 0) {
+ print_log(log_file, step, step * tau, prob);
+ }
+ }
+}
+
+int main(int argc, char* argv[]) {
+ if (argc < 4) {
+ std::cout << "Usage: " << argv[0] << "<heat conductivity> <nx> <nsteps_t> <dt>" << std::endl;
+ return 1;
+ }
+ double heat_cond = std::stod(argv[1]);
+ size_t nx = std::stoi(argv[2]);
+ size_t nsteps_t = std::stoi(argv[3]);
+ double dt = std::stod(argv[4]);
+
+ std::string log_name = "heat_eq_cg_log.a" +
+ std::to_string(heat_cond) +
+ ".grid" + std::to_string(nx) +
+ "x1.dt" + std::to_string(dt) +
+ ".txt";
+
+ std::ofstream log_file(log_name);
+ size_t log_freq = nsteps_t / 50;
+ if (log_freq == 0) {
+ log_freq = 1;
+ }
+
+ vec<double> us(nx+2), us_next(nx+2);
+ heat_eq_problem prob = {us, us_next, 1.0 / nx, heat_cond};
+
+ set_constant_u(prob.us, 1.0);
+ apply_dirichlet_bc(prob, 1.0, 0.5);
+
+ double t0 = omp_get_wtime();
+ dirichlet_implicit_solver(prob, 1.0, 0.5, 1.0, dt, nsteps_t, log_freq, log_file);
+ double t1 = omp_get_wtime();
+ std::cout << "Execution time: " << t1 - t0 << "\n";
+
+ return 0;
+}
diff --git a/heat_equation_explicit.cpp b/heat_equation_explicit.cpp
new file mode 100644
index 0000000..1043a07
--- /dev/null
+++ b/heat_equation_explicit.cpp
@@ -0,0 +1,156 @@
+#include "../vectors_and_matrices/array_types.hpp"
+#include "heat_equation.hpp"
+#include <iostream>
+#include <fstream>
+#include <cstdint>
+#include <cmath>
+#include <string>
+#include <omp.h>
+
+using intptr_t = std::intptr_t;
+using std::cin;
+using std::cout;
+using std::endl;
+using std::floor;
+
+void set_constant_u(vec<double> us, double u0)
+{
+ intptr_t n = us.length()-2;
+ for (intptr_t i=1; i<=n; i++)
+ {
+ us(i) = u0;
+ }
+}
+
+void apply_bc_explicit(heat_eq_problem& prob)
+{
+ boundary_cond left_bc = prob.left_bc, right_bc = prob.right_bc;
+
+ // Граничное условие имеет вид
+ // a u + b du/dx = z
+ // Аппроксимируем как
+ // a (u_0 + u_1) / 2 + b (u_1 - u_0) / h = z
+ vec<double> us = prob.us;
+ intptr_t n = us.length() - 2;
+ double h = prob.step_x;
+
+ double al = left_bc.coeff_u, bl = left_bc.coeff_du, zl = left_bc.rhp_val;
+ us(0) = (zl - (al / 2 + bl / h) * us(1)) / (al / 2 - bl / h);
+
+ double ar = right_bc.coeff_u, br = right_bc.coeff_du, zr = right_bc.rhp_val;
+ us(n+1) = (zr - (ar / 2 - br / h) * us(n)) / (ar / 2 + br / h);
+}
+
+void explicit_step(heat_eq_problem& prob, double tau)
+{
+ vec<double> us = prob.us;
+ vec<double> us_next = prob.us_next;
+ intptr_t n = us.length() - 2;
+ double h = prob.step_x;
+ double courant = tau * prob.conductivity / (h * h);
+
+ apply_bc_explicit(prob);
+ for (intptr_t i=1; i<=n; i++)
+ {
+ us_next(i) = us(i) + courant * (-2 * us(i) + us(i-1) + us(i+1));
+ }
+ // prob.us, prob.us_next = prob.us_next, prob.us
+ vec<double> tmp = us;
+ prob.us = us_next;
+ prob.us_next = tmp;
+}
+
+void explicit_solver(heat_eq_problem& prob,
+ double L, double tau, intptr_t ntimesteps,
+ intptr_t log_freq, std::ostream& log_file)
+{
+ // входные параметры:
+ // prob - начальная конфигурация
+ // L - длина отрезка, на котором решается уравнение
+ // tau - шаг по времени
+ // ntimesteps - сколько шагов по времени сделать
+ // log_freq - раз в столько шагов выводить информацию
+ // log_file - файловый дескриптор для вывода информации
+ intptr_t nx = prob.us.length() - 2;
+ double h = L / nx;
+ prob.step_x = h;
+ double courant = 2.0 * prob.conductivity * tau / (h * h);
+ if (courant > 1.0) {
+ std::cerr << "Warning: Courant number > 1\n";
+ log_file << "Warning: Courant number > 1\n";
+ }
+ apply_bc_explicit(prob);
+ print_log(log_file, 0, 0.0, prob);
+ // эволюция по времени
+ for (intptr_t step = 1; step <= ntimesteps; step++) {
+ explicit_step(prob, tau);
+ if (step % log_freq == 0) {
+ print_log(log_file, step, step * tau, prob);
+ }
+ }
+}
+
+int main(int argc, char* argv[]) {
+
+ std::string npoints_str, dx_str, dt_str, tmax_str;
+
+ double heat_cond = 1.0;
+ intptr_t nx;
+ double dx, dt, tmax;
+
+ cout << "Enter:\nNPOINTS: n\nSTEP_X: h\nSTEP_T: dt\nT_MAX: T" << endl;
+ cin >> npoints_str;
+ if (npoints_str != "NPOINTS:") {
+ cout << "ERROR: Cannot find NPOINTS record\n";
+ return 1;
+ }
+ cin >> nx;
+
+ cin >> dx_str;
+ if (dx_str != "STEP_X:") {
+ cout << "ERROR: Cannot find STEP_X record\n";
+ return 1;
+ }
+ cin >> dx;
+
+ cin >> dt_str;
+ if (dt_str != "STEP_T:") {
+ cout << "ERROR: Cannot find STEP_T record\n";
+ return 1;
+ }
+ cin >> dt;
+
+ cin >> tmax_str;
+ if (tmax_str != "T_MAX:") {
+ cout << "ERROR: Cannot find T_MAX record\n";
+ return 1;
+ }
+ cin >> tmax;
+
+ intptr_t nsteps_t = floor(tmax / dt);
+ std::string log_name = "heat_eq_expl_log.grid" + std::to_string(nx) + "pts.dx" + std::to_string(dx) + ".dt" + std::to_string(dt) + ".txt";
+
+ std::ofstream log_file(log_name);
+ intptr_t log_freq = nsteps_t / 5;
+ if (log_freq == 0) {
+ log_freq = 1;
+ }
+
+ vec<double> us(nx+2), us_next(nx+2);
+ heat_eq_problem prob = {us, us_next, dx, heat_cond, dirichlet_bc(1.0), dirichlet_bc(0.5)};
+
+ set_constant_u(prob.us, 1.0);
+
+ double t0 = omp_get_wtime();
+ explicit_solver(prob, nx * dx, dt, nsteps_t, log_freq, log_file);
+ double t1 = omp_get_wtime();
+
+ intptr_t nflops = (5 * nx + 18) * nsteps_t; // 5nx - основные итерации, 18 - применение ГУ
+
+ cout << "Execution time: " << t1 - t0 << "\n";
+ cout << "Steps per second: " << nsteps_t / (t1 - t0) << "\n";
+ cout << "GFLOPS: " << 1e-9 * nflops / (t1 - t0) << "\n";
+ cout << "Output written to " + log_name << std::endl;
+
+ return 0;
+}
diff --git a/heat_equation_explicit_omp.cpp b/heat_equation_explicit_omp.cpp
new file mode 100644
index 0000000..ed7d5ac
--- /dev/null
+++ b/heat_equation_explicit_omp.cpp
@@ -0,0 +1,157 @@
+#include "../vectors_and_matrices/array_types.hpp"
+#include "heat_equation.hpp"
+#include <iostream>
+#include <fstream>
+#include <cstdint>
+#include <cmath>
+#include <string>
+#include <omp.h>
+
+using intptr_t = std::intptr_t;
+using std::cin;
+using std::cout;
+using std::endl;
+using std::floor;
+
+void set_constant_u(vec<double> us, double u0)
+{
+ intptr_t n = us.length()-2;
+ for (intptr_t i=1; i<=n; i++)
+ {
+ us(i) = u0;
+ }
+}
+
+void apply_bc_explicit(heat_eq_problem& prob)
+{
+ boundary_cond left_bc = prob.left_bc, right_bc = prob.right_bc;
+
+ // Граничное условие имеет вид
+ // a u + b du/dx = z
+ // Аппроксимируем как
+ // a (u_0 + u_1) / 2 + b (u_1 - u_0) / h = z
+ vec<double> us = prob.us;
+ intptr_t n = us.length() - 2;
+ double h = prob.step_x;
+
+ double al = left_bc.coeff_u, bl = left_bc.coeff_du, zl = left_bc.rhp_val;
+ us(0) = (zl - (al / 2 + bl / h) * us(1)) / (al / 2 - bl / h);
+
+ double ar = right_bc.coeff_u, br = right_bc.coeff_du, zr = right_bc.rhp_val;
+ us(n+1) = (zr - (ar / 2 - br / h) * us(n)) / (ar / 2 + br / h);
+}
+
+void explicit_step(heat_eq_problem& prob, double tau)
+{
+ vec<double> us = prob.us;
+ vec<double> us_next = prob.us_next;
+ intptr_t n = us.length() - 2;
+ double h = prob.step_x;
+ double courant = tau * prob.conductivity / (h * h);
+
+ apply_bc_explicit(prob);
+ #pragma omp parallel for
+ for (intptr_t i=1; i<=n; i++)
+ {
+ us_next(i) = us(i) + courant * (-2 * us(i) + us(i-1) + us(i+1));
+ }
+ // prob.us, prob.us_next = prob.us_next, prob.us
+ vec<double> tmp = us;
+ prob.us = us_next;
+ prob.us_next = tmp;
+}
+
+void explicit_solver(heat_eq_problem& prob,
+ double L, double tau, intptr_t ntimesteps,
+ intptr_t log_freq, std::ostream& log_file)
+{
+ // входные параметры:
+ // prob - начальная конфигурация
+ // L - длина отрезка, на котором решается уравнение
+ // tau - шаг по времени
+ // ntimesteps - сколько шагов по времени сделать
+ // log_freq - раз в столько шагов выводить информацию
+ // log_file - файловый дескриптор для вывода информации
+ intptr_t nx = prob.us.length() - 2;
+ double h = L / nx;
+ prob.step_x = h;
+ double courant = 2.0 * prob.conductivity * tau / (h * h);
+ if (courant > 1.0) {
+ std::cerr << "Warning: Courant number > 1\n";
+ log_file << "Warning: Courant number > 1\n";
+ }
+ apply_bc_explicit(prob);
+ print_log(log_file, 0, 0.0, prob);
+ // эволюция по времени
+ for (intptr_t step = 1; step <= ntimesteps; step++) {
+ explicit_step(prob, tau);
+ if (step % log_freq == 0) {
+ print_log(log_file, step, step * tau, prob);
+ }
+ }
+}
+
+int main(int argc, char* argv[]) {
+
+ std::string npoints_str, dx_str, dt_str, tmax_str;
+
+ double heat_cond = 1.0;
+ intptr_t nx;
+ double dx, dt, tmax;
+
+ cout << "Enter:\nNPOINTS: n\nSTEP_X: h\nSTEP_T: dt\nT_MAX: T" << endl;
+ cin >> npoints_str;
+ if (npoints_str != "NPOINTS:") {
+ cout << "ERROR: Cannot find NPOINTS record\n";
+ return 1;
+ }
+ cin >> nx;
+
+ cin >> dx_str;
+ if (dx_str != "STEP_X:") {
+ cout << "ERROR: Cannot find STEP_X record\n";
+ return 1;
+ }
+ cin >> dx;
+
+ cin >> dt_str;
+ if (dt_str != "STEP_T:") {
+ cout << "ERROR: Cannot find STEP_T record\n";
+ return 1;
+ }
+ cin >> dt;
+
+ cin >> tmax_str;
+ if (tmax_str != "T_MAX:") {
+ cout << "ERROR: Cannot find T_MAX record\n";
+ return 1;
+ }
+ cin >> tmax;
+
+ intptr_t nsteps_t = floor(tmax / dt);
+ std::string log_name = "heat_eq_expl_log.grid" + std::to_string(nx) + "pts.dx" + std::to_string(dx) + ".dt" + std::to_string(dt) + ".txt";
+
+ std::ofstream log_file(log_name);
+ intptr_t log_freq = nsteps_t / 5;
+ if (log_freq == 0) {
+ log_freq = 1;
+ }
+
+ vec<double> us(nx+2), us_next(nx+2);
+ heat_eq_problem prob = {us, us_next, dx, heat_cond, dirichlet_bc(1.0), dirichlet_bc(0.5)};
+
+ set_constant_u(prob.us, 1.0);
+
+ double t0 = omp_get_wtime();
+ explicit_solver(prob, nx * dx, dt, nsteps_t, log_freq, log_file);
+ double t1 = omp_get_wtime();
+
+ intptr_t nflops = (5 * nx + 18) * nsteps_t; // 5nx - основные итерации, 18 - применение ГУ
+
+ cout << "Execution time: " << t1 - t0 << "\n";
+ cout << "Steps per second: " << nsteps_t / (t1 - t0) << "\n";
+ cout << "GFLOPS: " << 1e-9 * nflops / (t1 - t0) << "\n";
+ cout << "Output written to " + log_name << std::endl;
+
+ return 0;
+}
diff --git a/heat_equation_implicit.cpp b/heat_equation_implicit.cpp
new file mode 100644
index 0000000..bc931e3
--- /dev/null
+++ b/heat_equation_implicit.cpp
@@ -0,0 +1,188 @@
+#include "../vectors_and_matrices/array_types.hpp"
+#include "heat_equation.hpp"
+#include <iostream>
+#include <fstream>
+#include <cstdint>
+#include <cmath>
+#include <string>
+#include <omp.h>
+
+using intptr_t = std::intptr_t;
+using std::cin;
+using std::cout;
+using std::endl;
+using std::floor;
+
+void set_constant_u(vec<double> us, double u0)
+{
+ intptr_t n = us.length()-2;
+ for (intptr_t i=1; i<=n; i++)
+ {
+ us(i) = u0;
+ }
+}
+
+void solve_tridiag(vec<double> x,
+ vec<double> lower, vec<double> diag, vec<double> upper,
+ vec<double> rhp
+ )
+{
+ // входные параметры:
+ // lower, diag, upper — массивы коэффициентов l, d, u
+ // rhp — массив с правой частью (right-hand part)
+ // x — массив, в который будет записываться решение
+ // courant - a * tau / h^2
+
+ // Замечание: если rhp и x — это один и тот же массив, то программа корректно отработает, только правая часть будет перезаписана решением
+ // РїСЂСЏРјРѕР№ С…РѕРґ
+ double factor;
+ intptr_t n = rhp.length()-2;
+ x(0) = rhp(0);
+ for (intptr_t i = 1; i <= n+1; i++)
+ {
+ factor = lower(i) / diag(i-1);
+ diag(i) -= factor * upper(i-1);
+ x(i) = rhp(i) - factor * x(i-1);
+ }
+ x(n+1) /= diag(n+1);
+ // обратный ход
+ for (intptr_t i = n+1; i > 0; i--)
+ {
+ x(i-1) = (x(i-1) - upper(i-1) * x(i)) / diag(i-1);
+ }
+}
+
+void apply_bc_implicit(heat_eq_problem& prob,
+ vec<double> lower, vec<double> diag, vec<double> upper)
+{
+ boundary_cond left_bc = prob.left_bc, right_bc = prob.right_bc;
+
+ // Граничное условие имеет вид
+ // a u + b du/dx = z
+ // Аппроксимируем как
+ // a (u_0 + u_1) / 2 + b (u_1 - u_0) / h = z
+ vec<double> us = prob.us;
+ intptr_t n = us.length() - 2;
+ double h = prob.step_x;
+
+ double al = left_bc.coeff_u, bl = left_bc.coeff_du, zl = left_bc.rhp_val;
+ diag(0) = (al / 2 - bl / h);
+ upper(0) = (al / 2 + bl / h);
+ us(0) = zl;
+
+ double ar = right_bc.coeff_u, br = right_bc.coeff_du, zr = right_bc.rhp_val;
+ diag(n+1) = (ar / 2 + br / h);
+ lower(n+1) = (ar / 2 - br / h);
+ us(n+1) = zr;
+}
+
+void implicit_step(heat_eq_problem& prob, double tau,
+ vec<double> lower, vec<double> diag, vec<double> upper)
+{
+ vec<double> us = prob.us;
+ intptr_t n = us.length() - 2;
+ double h = prob.step_x;
+ double courant = tau * prob.conductivity / (h * h);
+
+ set_constant_u(diag, 1.0 + 2 * courant);
+
+ apply_bc_implicit(prob, lower, diag, upper);
+ solve_tridiag(us, lower, diag, upper, us);
+}
+
+void implicit_solver(heat_eq_problem& prob,
+ double L, double tau, intptr_t ntimesteps,
+ intptr_t log_freq, std::ostream& log_file)
+{
+ // входные параметры:
+ // prob - начальная конфигурация
+ // L - длина отрезка, на котором решается уравнение
+ // tau - шаг по времени
+ // ntimesteps - сколько шагов по времени сделать
+ // log_freq - раз в столько шагов выводить информацию
+ // log_file - файловый дескриптор для вывода информации
+ intptr_t nx = prob.us.length() - 2;
+ double h = L / nx;
+ prob.step_x = h;
+ double courant = prob.conductivity * tau / (h * h);
+ if (courant > 0.5) {
+ std::cerr << "Warning: Courant number > 1\n";
+ }
+ print_log(log_file, 0, 0.0, prob);
+ vec<double> lower(nx+2);
+ vec<double> diag(nx+2);
+ vec<double> upper(nx+2);
+ set_constant_u(lower, -courant);
+ set_constant_u(upper, -courant);
+ // эволюция по времени
+ for (intptr_t step = 1; step <= ntimesteps; step++) {
+ implicit_step(prob, tau, lower, diag, upper);
+ if (step % log_freq == 0) {
+ print_log(log_file, step, step * tau, prob);
+ }
+ }
+}
+
+int main(int argc, char* argv[]) {
+ std::string npoints_str, dx_str, dt_str, tmax_str;
+
+ double heat_cond = 1.0;
+ intptr_t nx;
+ double dx, dt, tmax;
+
+ cout << "Enter:\nNPOINTS: n\nSTEP_X: h\nSTEP_T: dt\nT_MAX: T" << endl;
+ cin >> npoints_str;
+ if (npoints_str != "NPOINTS:") {
+ cout << "ERROR: Cannot find NPOINTS record\n";
+ return 1;
+ }
+ cin >> nx;
+
+ cin >> dx_str;
+ if (dx_str != "STEP_X:") {
+ cout << "ERROR: Cannot find STEP_X record\n";
+ return 1;
+ }
+ cin >> dx;
+
+ cin >> dt_str;
+ if (dt_str != "STEP_T:") {
+ cout << "ERROR: Cannot find STEP_T record\n";
+ return 1;
+ }
+ cin >> dt;
+
+ cin >> tmax_str;
+ if (tmax_str != "T_MAX:") {
+ cout << "ERROR: Cannot find T_MAX record\n";
+ return 1;
+ }
+ cin >> tmax;
+
+ intptr_t nsteps_t = floor(tmax / dt);
+ std::string log_name = "heat_eq_impl_log.grid" + std::to_string(nx) + "pts.dx" + std::to_string(dx) + ".dt" + std::to_string(dt) + ".txt";
+
+ std::ofstream log_file(log_name);
+ std::intptr_t log_freq = nsteps_t / 5;
+ if (log_freq == 0) {
+ log_freq = 1;
+ }
+
+ vec<double> us(nx+2), us_next(nx+2);
+ heat_eq_problem prob = {us, us_next, dx, heat_cond, dirichlet_bc(1.0), dirichlet_bc(0.5)};
+
+ set_constant_u(prob.us, 1.0);
+
+ double t0 = omp_get_wtime();
+ implicit_solver(prob, nx * dx, dt, nsteps_t, log_freq, log_file);
+ double t1 = omp_get_wtime();
+
+ intptr_t nflops = (8 * nx + 12) * nsteps_t; // 8nx - прогонка, 12 - применение ГУ
+
+ cout << "Execution time: " << t1 - t0 << "\n";
+ cout << "Steps per second: " << nsteps_t / (t1 - t0) << "\n";
+ cout << "GFLOPS: " << 1e-9 * nflops / (t1 - t0) << "\n";
+ cout << "Output written to " + log_name << std::endl;
+
+ return 0;
+}
--
GitLab