diffusion.h
I'm done with the backbone of it, I haven't had feedback on it.
This commit is contained in:
Executable
BIN
Binary file not shown.
@@ -0,0 +1,36 @@
|
||||
#pragma once
|
||||
|
||||
|
||||
namespace core{
|
||||
|
||||
enum class GridKind { Uniform, NonUniform };
|
||||
enum class FDKind { Central, Forward, Backward };
|
||||
enum class BCKind { Dirichlet, Neumann /*, Robin*/ };
|
||||
enum class SolverKind { LU, Inverse /*, CG*/ };
|
||||
|
||||
template <typename T>
|
||||
struct BC {
|
||||
FDKind fd{FDKind::Forward};
|
||||
BCKind kind{BCKind::Dirichlet};
|
||||
T value{T(0)};
|
||||
};
|
||||
|
||||
|
||||
// Global default config holder
|
||||
template <typename T>
|
||||
struct Configs {
|
||||
GridKind grid{GridKind::Uniform};
|
||||
FDKind fd{FDKind::Central};
|
||||
BC<T> left{FDKind::Forward, BCKind::Dirichlet, T(0) };
|
||||
BC<T> right{FDKind::Backward, BCKind::Dirichlet, T(0) };
|
||||
SolverKind solver{SolverKind::LU};
|
||||
|
||||
static Configs& defaults() {
|
||||
static Configs g{}; // process-wide defaults
|
||||
return g;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
|
||||
} // namespace core
|
||||
@@ -1,6 +1,6 @@
|
||||
|
||||
|
||||
#pragma once
|
||||
|
||||
#include <vector>
|
||||
#include <omp.h>
|
||||
|
||||
|
||||
|
||||
@@ -0,0 +1,129 @@
|
||||
#pragma once
|
||||
|
||||
#include "modules/mesh/mesh1d.h"
|
||||
#include "core/global_config.h"
|
||||
#include "utils/matrix.h"
|
||||
#include "utils/vector.h"
|
||||
|
||||
|
||||
|
||||
|
||||
namespace fluids {
|
||||
|
||||
template <typename T>
|
||||
struct Diffusion1D{
|
||||
const core::Configs<T>& cfg;
|
||||
const mesh::Mesh1D<T>& mesh;
|
||||
T Gamma{1};
|
||||
|
||||
|
||||
|
||||
// Constructor
|
||||
Diffusion1D(const core::Configs<T>& configs, const mesh::Mesh1D<T>& Mesh, T Gamma_const=T(1)): cfg(configs), mesh(Mesh), Gamma(Gamma_const) {}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
void assemble(utils::Matrix<T>& A, utils::Vector<T>& b, utils::Vector<T>& s){
|
||||
|
||||
uint64_t N = mesh.center_idx + 1;
|
||||
|
||||
if (N < 3){
|
||||
throw std::runtime_error("Diffusion1D: need N>=3");
|
||||
}
|
||||
if (A.rows() != N || A.cols() != N){
|
||||
A = utils::Matrix<T>(N, N, T(0));
|
||||
}
|
||||
if (b.size() != N){
|
||||
b = utils::Vector<T>(N, T(0));
|
||||
}
|
||||
|
||||
|
||||
if (cfg.grid == core::GridKind::Uniform){
|
||||
// Core of A
|
||||
if (cfg.fd == core::FDKind::Central){
|
||||
uniform_central_finite_diffrence_2_order(A,b,s);
|
||||
}
|
||||
|
||||
|
||||
|
||||
// Left BC of A
|
||||
if (cfg.left.kind == core::BCKind::Dirichlet){
|
||||
BC_uniform_backward_finite_diffrence_2_order_Dirichlet(A,b,s);
|
||||
}else if (cfg.left.kind == core::BCKind::Neumann){
|
||||
BC_uniform_backward_finite_diffrence_2_order_Neumann(A,b,s);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void uniform_central_finite_diffrence_2_order(utils::Matrix<T>& A, utils::Vector<T>& b, utils::Vector<T>& s){
|
||||
T xm, xc, xp;
|
||||
|
||||
for (uint64_t i = 1; i < mesh.center_idx; ++i){
|
||||
xm = mesh.center(i-1);
|
||||
xc = mesh.center(i);
|
||||
xp = mesh.center(i+1);
|
||||
|
||||
A(i, i-1) = Gamma/(xc - xm);
|
||||
A(i, i) = -((Gamma/(xp - xc)) + (Gamma/(xc - xm)));
|
||||
A(i, i+1) = Gamma/(xp - xc);
|
||||
|
||||
b[i] = -s[i]*mesh.dx(i);
|
||||
}
|
||||
}
|
||||
|
||||
void BC_uniform_backward_finite_diffrence_2_order_Dirichlet(utils::Matrix<T>& A, utils::Vector<T>& b, utils::Vector<T>& s){
|
||||
T xm;
|
||||
T xw = mesh.vertice(0);
|
||||
T xc = mesh.center(0);
|
||||
T xp = mesh.center(1);
|
||||
T xe;
|
||||
uint64_t N = mesh.center_idx;
|
||||
|
||||
A(0, 0) = -((Gamma/(xp - xc)) + (Gamma/(xc - xw)));
|
||||
A(0, 1) = Gamma/(xp - xc);
|
||||
|
||||
b[0] = -s[0]*mesh.dx(0) - Gamma*(cfg.left.value/(xc - xw));
|
||||
|
||||
xm = mesh.center(N-1);
|
||||
xc = mesh.center(N);
|
||||
xe = mesh.vertice(N+1);
|
||||
|
||||
A(N, N-1) = Gamma/(xc - xm);
|
||||
A(N, N) = -((Gamma/(xe - xc)) + (Gamma/(xc - xm)));
|
||||
|
||||
b[N] = -s[N]*mesh.dx(N) - Gamma*(cfg.right.value/(xe - xc));
|
||||
A.print();
|
||||
b.print();
|
||||
}
|
||||
|
||||
void BC_uniform_backward_finite_diffrence_2_order_Neumann(utils::Matrix<T>& A, utils::Vector<T>& b, utils::Vector<T>& s){
|
||||
T xm;
|
||||
T xc = mesh.center(0);
|
||||
T xp = mesh.center(1);
|
||||
|
||||
uint64_t N = mesh.center_idx;
|
||||
|
||||
A(0, 0) = -Gamma/(xp - xc);
|
||||
A(0, 1) = Gamma/(xp - xc);
|
||||
|
||||
b[0] = -s[0]*mesh.dx(0) - (Gamma*cfg.left.value);
|
||||
|
||||
xm = mesh.center(N-1);
|
||||
xc = mesh.center(N);
|
||||
|
||||
A(N, N-1) = Gamma/(xc - xm);
|
||||
A(N, N) = -Gamma/(xc - xm);
|
||||
|
||||
b[N] = -s[N]*mesh.dx(N) - Gamma*(cfg.right.value);
|
||||
|
||||
A.print();
|
||||
b.print();
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
} // end namespace fluids
|
||||
@@ -0,0 +1,5 @@
|
||||
#pragma once
|
||||
|
||||
#include "modules/fluids/diffusion1d.h"
|
||||
|
||||
|
||||
@@ -0,0 +1,4 @@
|
||||
#pragma once
|
||||
|
||||
#include "modules/mesh/mesh1d.h"
|
||||
|
||||
@@ -2,20 +2,18 @@
|
||||
|
||||
#include "utils/vector.h"
|
||||
|
||||
|
||||
namespace fvm {
|
||||
namespace mesh {
|
||||
|
||||
template <typename T>
|
||||
|
||||
struct Grid1D{
|
||||
struct Mesh1D{
|
||||
uint64_t center_idx; // max cell index
|
||||
uint64_t vertices_idx; // max vertice index
|
||||
utils::Vector<T> centers; // size N (unknowns at cell centers)
|
||||
utils::Vector<T> vertices; // size N+1 (face positions)
|
||||
|
||||
Grid1D() = default;
|
||||
Mesh1D() = default;
|
||||
|
||||
explicit Grid1D(const utils::Vector<T>& midpoints){
|
||||
explicit Mesh1D(const utils::Vector<T>& midpoints){
|
||||
centers = midpoints;
|
||||
center_idx = centers.size()-1;
|
||||
vertices_idx = centers.size();
|
||||
@@ -33,11 +31,12 @@ namespace fvm {
|
||||
|
||||
T dx(uint64_t i) const { check(i); return vertices[i+1] - vertices[i]; }
|
||||
T center(uint64_t i) const { check(i); return centers[i]; }
|
||||
T vertice(uint64_t i) const {; return vertices[i]; }
|
||||
|
||||
void check(uint64_t i) const {
|
||||
if (i > center_idx) throw std::runtime_error("Grid1D: cell index out of range");
|
||||
if (i > center_idx) throw std::runtime_error("Mesh1D: cell index out of range");
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
} // end namespace mesh
|
||||
@@ -10,10 +10,6 @@
|
||||
#include "./numerics/min.h"
|
||||
#include "./numerics/max.h"
|
||||
#include "./numerics/abs.h"
|
||||
#include "./numerics/interpolation1d_base.h" // base
|
||||
#include "./numerics/interpolation1d/interpolation1d_linear.h" // derived
|
||||
#include "./numerics/interpolation1d/interpolation1d_polynomial.h" // derived
|
||||
#include "./numerics/interpolation1d/interpolation1d_cubic_spline.h" // derived
|
||||
#include "./numerics/interpolation1d/interpolation1d_rational.h" // derived
|
||||
#include "./numerics/interpolation1d/interpolation1d_barycentric.h" // derived
|
||||
#include "./numerics/interpolation1d.h" // base
|
||||
|
||||
|
||||
|
||||
@@ -0,0 +1,4 @@
|
||||
// "./utils/generators.h"
|
||||
#pragma once
|
||||
|
||||
#include "./utils/generators/linspace.h"
|
||||
@@ -0,0 +1,37 @@
|
||||
#pragma once
|
||||
|
||||
#include "utils/vector.h"
|
||||
|
||||
|
||||
namespace utils{
|
||||
|
||||
template <typename T>
|
||||
void inplace_linspace(utils::Vector<T>& a, T start, T stop, bool endpoint=true){
|
||||
|
||||
uint64_t N = a.size();
|
||||
T step;
|
||||
|
||||
if (endpoint){
|
||||
step = (stop - start) / static_cast<T>(N - 1);
|
||||
}else{
|
||||
step = (stop - start) / static_cast<T>(N);
|
||||
}
|
||||
|
||||
for (uint64_t i = 0; i < N; ++i){
|
||||
a[i] = start + (step*static_cast<T>(i));
|
||||
}
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
utils::Vector<T> linspace(T start, T stop, uint64_t N, bool endpoint=true){
|
||||
|
||||
utils::Vector<T> a(N);
|
||||
|
||||
inplace_linspace(a, start, stop, endpoint);
|
||||
|
||||
return a;
|
||||
}
|
||||
|
||||
|
||||
|
||||
} // end namespace utils
|
||||
@@ -3,3 +3,4 @@
|
||||
|
||||
#include "./utils/vector.h"
|
||||
#include "./utils/matrix.h"
|
||||
#include "./utils/generators.h"
|
||||
+56
@@ -0,0 +1,56 @@
|
||||
obj/main.o: src/main.cpp include/./core/omp_config.h \
|
||||
include/./utils/utils.h include/./utils/vector.h \
|
||||
include/./utils/matrix.h include/./utils/generators.h \
|
||||
include/./utils/generators/linspace.h include/utils/vector.h \
|
||||
include/./numerics/numerics.h include/./numerics/initializers/eye.h \
|
||||
include/./numerics/matequal.h include/./numerics/abs.h \
|
||||
include/./numerics/transpose.h include/./numerics/inverse.h \
|
||||
include/./numerics/inverse/inverse_gauss_jordan.h \
|
||||
include/./numerics/inverse/inverse_lu.h include/./decomp/lu.h \
|
||||
include/./numerics/matmul.h include/./numerics/matvec.h \
|
||||
include/./numerics/min.h include/./numerics/max.h \
|
||||
include/./numerics/interpolation1d.h \
|
||||
include/./numerics/interpolation1d/interpolation1d_barycentric.h \
|
||||
include/./numerics/interpolation1d/interpolation1d_base.h \
|
||||
include/./numerics/interpolation1d/interpolation1d_cubic_spline.h \
|
||||
include/./numerics/interpolation1d/interpolation1d_linear.h \
|
||||
include/./numerics/interpolation1d/interpolation1d_polynomial.h \
|
||||
include/./numerics/interpolation1d/interpolation1d_rational.h \
|
||||
include/./decomp/decomp.h include/./modules/mesh/mesh.h \
|
||||
include/modules/mesh/mesh1d.h include/modules/fluids/fluids.h \
|
||||
include/modules/fluids/diffusion1d.h include/core/global_config.h \
|
||||
include/utils/matrix.h
|
||||
include/./core/omp_config.h:
|
||||
include/./utils/utils.h:
|
||||
include/./utils/vector.h:
|
||||
include/./utils/matrix.h:
|
||||
include/./utils/generators.h:
|
||||
include/./utils/generators/linspace.h:
|
||||
include/utils/vector.h:
|
||||
include/./numerics/numerics.h:
|
||||
include/./numerics/initializers/eye.h:
|
||||
include/./numerics/matequal.h:
|
||||
include/./numerics/abs.h:
|
||||
include/./numerics/transpose.h:
|
||||
include/./numerics/inverse.h:
|
||||
include/./numerics/inverse/inverse_gauss_jordan.h:
|
||||
include/./numerics/inverse/inverse_lu.h:
|
||||
include/./decomp/lu.h:
|
||||
include/./numerics/matmul.h:
|
||||
include/./numerics/matvec.h:
|
||||
include/./numerics/min.h:
|
||||
include/./numerics/max.h:
|
||||
include/./numerics/interpolation1d.h:
|
||||
include/./numerics/interpolation1d/interpolation1d_barycentric.h:
|
||||
include/./numerics/interpolation1d/interpolation1d_base.h:
|
||||
include/./numerics/interpolation1d/interpolation1d_cubic_spline.h:
|
||||
include/./numerics/interpolation1d/interpolation1d_linear.h:
|
||||
include/./numerics/interpolation1d/interpolation1d_polynomial.h:
|
||||
include/./numerics/interpolation1d/interpolation1d_rational.h:
|
||||
include/./decomp/decomp.h:
|
||||
include/./modules/mesh/mesh.h:
|
||||
include/modules/mesh/mesh1d.h:
|
||||
include/modules/fluids/fluids.h:
|
||||
include/modules/fluids/diffusion1d.h:
|
||||
include/core/global_config.h:
|
||||
include/utils/matrix.h:
|
||||
BIN
Binary file not shown.
+23
-123
@@ -1,146 +1,46 @@
|
||||
#include "./core/omp_config.h"
|
||||
|
||||
#include "./utils/utils.h"
|
||||
#include "./numerics/numerics.h"
|
||||
#include "./decomp/decomp.h"
|
||||
#include "./core/omp_config.h"
|
||||
|
||||
#include "./modules/grid1d.h"
|
||||
#include "./modules/finitedifference1d.h"
|
||||
|
||||
|
||||
#include <iostream>
|
||||
#include <stdexcept>
|
||||
#include "./modules/mesh/mesh.h"
|
||||
#include "modules/fluids/fluids.h"
|
||||
|
||||
|
||||
#include <chrono>
|
||||
|
||||
//#include "./numerics/interpolation/interpolation_linear.h"
|
||||
//#include <iostream>
|
||||
//#include <stdexcept>
|
||||
//#include <chrono>
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
int main(int argc, char const *argv[])
|
||||
{
|
||||
utils::Vd a = utils::linspace<double>(1, 10, 10, true);
|
||||
a.print();
|
||||
mesh::Mesh1D<double> mesh(a);
|
||||
mesh.generate_vertices(0.5, 10.5);
|
||||
double Gamma = 1.0;
|
||||
|
||||
|
||||
utils::Md A;
|
||||
|
||||
/*
|
||||
int hw = omp_get_max_active_levels();
|
||||
if (hw <= 1) return 0;
|
||||
|
||||
const uint64_t m=512, k=512, p=512; // ~134M MACs; adjust if needed
|
||||
utils::Md A(m,k,1), B(k,p,1), C(k,p,1);
|
||||
|
||||
omp_set_max_active_levels(1);
|
||||
|
||||
auto t0 = std::chrono::high_resolution_clock::now();
|
||||
for (int i = 0; i < m*k*p; ++i){
|
||||
A==B
|
||||
}
|
||||
double t1 = std::chrono::duration<double>(std::chrono::high_resolution_clock::now() - t0).count();
|
||||
utils::Vd b, s(10,1);
|
||||
|
||||
|
||||
omp_set_max_active_levels(2);
|
||||
auto t0 = std::chrono::high_resolution_clock::now();
|
||||
for (int i = 0; i < m*k*p; ++i){
|
||||
A==B
|
||||
}
|
||||
double t1 = std::chrono::duration<double>(std::chrono::high_resolution_clock::now() - t0).count();
|
||||
|
||||
omp_set_num_threads(prev);
|
||||
|
||||
// Must not be notably slower with many threads
|
||||
CHECK(tN <= t1 * 1.05, "rows_omp: multi-thread slower than single-thread");
|
||||
core::Configs<double>& cfg = core::Configs<double>::defaults();
|
||||
cfg.grid = core::GridKind::Uniform;
|
||||
cfg.left = {core::FDKind::Forward, core::BCKind::Neumann, 0.0};
|
||||
cfg.right = {core::FDKind::Backward, core::BCKind::Neumann, 0.0};
|
||||
cfg.solver = core::SolverKind::LU;
|
||||
|
||||
|
||||
|
||||
utils::Md A(5,5, 1);
|
||||
utils::Md B(5,5, 1);
|
||||
utils::Md C(5,5, 2);
|
||||
|
||||
bool result1 = (A==B);
|
||||
bool result2 = (A==C);
|
||||
|
||||
omp_set_max_active_levels(1):
|
||||
|
||||
for (int i = 0; i < 100; ++i){
|
||||
(A==B)
|
||||
}
|
||||
|
||||
omp_set_max_active_levels(2):
|
||||
|
||||
for (int i = 0; i < 100; ++i){
|
||||
(A==B)
|
||||
}
|
||||
|
||||
std::cout << result1 << std::endl;
|
||||
std::cout << result2 << std::endl;
|
||||
fluids::Diffusion1D<double> diffusion(cfg, mesh, Gamma);
|
||||
diffusion.assemble(A, b, s);
|
||||
|
||||
|
||||
|
||||
|
||||
*/
|
||||
|
||||
|
||||
/*
|
||||
|
||||
utils::Vector<double> x(100, 0), y(100,0);
|
||||
for (uint64_t i = 0; i < 100; ++i){
|
||||
x[i] = i+1;
|
||||
y[i] = 1 + i*0.1;
|
||||
}
|
||||
//y[9] = 1.5;
|
||||
|
||||
x.print();
|
||||
y.print();
|
||||
|
||||
double p = 5.5;
|
||||
|
||||
|
||||
|
||||
numerics::interp_linear<double> linear(x,y);
|
||||
numerics::interp_polynomial<double> polynomial(x,y, 3);
|
||||
numerics::interp_cubic_spline<double> cubic_spline(x,y);
|
||||
numerics::interp_rational<double> rational(x,y,2);
|
||||
numerics::interp_barycentric<double> barycentric(x,y, 2);
|
||||
|
||||
std::cout << "=== interpolate: " << p << " ===" << std::endl;
|
||||
|
||||
std::cout << linear.interp(p) << std::endl;
|
||||
std::cout << linear.interp(p) << std::endl;
|
||||
std::cout << polynomial.interp(p) << std::endl; // error = polynomial.dy
|
||||
std::cout << cubic_spline.interp(p) << std::endl;
|
||||
std::cout << rational.interp(p) << std::endl;
|
||||
std::cout << barycentric.interp(p) << std::endl;
|
||||
|
||||
p += 0.01;
|
||||
|
||||
std::cout << "=== interpolate: " << p << " ===" << std::endl;
|
||||
|
||||
std::cout << linear.interp(p) << std::endl;
|
||||
std::cout << polynomial.interp(p) << std::endl;
|
||||
std::cout << cubic_spline.interp(p) << std::endl;
|
||||
std::cout << rational.interp(p) << std::endl;
|
||||
std::cout << barycentric.interp(p) << std::endl;
|
||||
|
||||
p += 0.01;
|
||||
|
||||
std::cout << "=== interpolate: " << p << " ===" << std::endl;
|
||||
|
||||
std::cout << linear.interp(p) << std::endl;
|
||||
std::cout << polynomial.interp(p) << std::endl;
|
||||
std::cout << cubic_spline.interp(p) << std::endl;
|
||||
std::cout << rational.interp(p) << std::endl;
|
||||
std::cout << barycentric.interp(p) << std::endl;
|
||||
|
||||
p += 50.01;
|
||||
|
||||
std::cout << "=== interpolate: " << p << " ===" << std::endl;
|
||||
|
||||
std::cout << linear.interp(p) << std::endl;
|
||||
std::cout << polynomial.interp(p) << std::endl;
|
||||
std::cout << cubic_spline.interp(p) << std::endl;
|
||||
std::cout << rational.interp(p) << std::endl;
|
||||
std::cout << barycentric.interp(p) << std::endl;
|
||||
|
||||
*/
|
||||
return 0;
|
||||
}
|
||||
Reference in New Issue
Block a user