Skip to content

Instantly share code, notes, and snippets.

@uranix

uranix/weno.cpp Secret

Created Sep 8, 2015
Embed
What would you like to do?
#include <iostream>
#include <cmath>
#include <cstdlib>
#include <vector>
template<typename T>
struct wrapped {
std::vector<T> v;
wrapped(std::size_t n, const T &e = T()) : v(n, e) { }
std::size_t size() const { return v.size(); }
std::ptrdiff_t idx(std::ptrdiff_t i) const {
const std::size_t s = size();
if (i < 0)
return i + s;
if (i >= s)
return i - s;
return i;
}
const T &operator[](std::ptrdiff_t i) const { return v[idx(i)]; }
T &operator[](std::ptrdiff_t i) { return v[idx(i)]; }
};
struct poly {
double a, b, c;
poly() { a = b = c = 0; }
poly(double a, double b, double c) : a(a), b(b), c(c) { }
double operator()(double x) const {
return a * x * x + b * x + c;
}
poly operator*(double v) const {
return poly(a * v, b * v, c * v);
}
poly operator+(const poly &o) const {
return poly(a + o.a, b + o.b, c + o.c);
}
poly &operator+=(const poly &o) {
a += o.a;
b += o.b;
c += o.c;
return *this;
}
static poly interpolate(double ul, double uc, double ur) {
double a = 0.5 * (ul - 2 * uc + ur);
double b = 0.5 * (ur - ul);
return poly(a, b, uc);
}
// Return q(x) = p(x - s)
poly shift(double s) {
return poly(a, b - 2 * a * s, (*this)(-s));
}
double smoothness() const {
return 2 * a + a * a / 3 + b * b;
}
};
struct poly5 {
double a, b, c, d, e;
poly5() {
a = b = c = d = e = 0;
}
void add(const poly &g, const poly &p) {
a += g.a * p.a;
b += g.a * p.b + g.b * p.a;
c += g.a * p.c + g.b * p.b + g.c * p.a;
d += g.b * p.c + g.c * p.b;
e += g.c * p.c;
}
double operator()(double x) const {
return e + x * (d + x * (c + x * (b + x * a)));
}
};
struct rat {
poly5 num;
poly denom;
rat(const poly5 &num = poly5(), const poly &denom = poly()) : num(num), denom(denom) { }
double operator()(const double x) {
return num(x) / denom(x);
}
};
template<typename Cont>
rat weno(const Cont &u, std::ptrdiff_t i) {
const poly pl = poly::interpolate(u[i-2], u[i-1], u[i]).shift(-1);
const poly pc = poly::interpolate(u[i-1], u[i], u[i+1]);
const poly pr = poly::interpolate(u[i], u[i+1], u[i+2]).shift(1);
const poly gl(1. / 12, -1. / 4, 1. / 6);
const poly gc(-1. / 6, 0, 2. / 3);
const poly gr(1. / 12, 1. / 4, 1. / 6);
poly5 num;
poly denom;
const double eps = 1e-6;
double al = std::pow(pl.smoothness() + eps, -2);
double ac = std::pow(pc.smoothness() + eps, -2);
double ar = std::pow(pr.smoothness() + eps, -2);
// al = ac = ar = 1;
num.add(gl * al, pl);
num.add(gc * ac, pc);
num.add(gr * ar, pr);
denom += gl * al;
denom += gc * ac;
denom += gr * ar;
return rat(num, denom);
}
double u0(double x) {
while (x < 0)
x += 1;
while (x > 1)
x -= 1;
return 0.001e-100 * sin(156 * x) + ((x > .2) && (x < .4) ? 1 : 0);
// return sin(2 * 3.14159265358929 * x);
}
int main() {
const size_t N = 100;
wrapped<double> u(N);
wrapped<rat> p(u.size());
const double h = 1. / N;
for (size_t i = 0; i < N; i++) {
double x = (i + .5) * h;
u[i] = u0(x);
}
const double C = .15;
const double dt = C * h;
double t = 0;
for (int steps = 0; steps < 20 * N; steps ++) {
for (size_t i = 0; i < N; i++)
p[i] = weno(u, i);
for (size_t i = 0; i < N; i++)
u[i] += C * (p[i-1](.5) - p[i](.5));
t += dt;
}
for (size_t i = 0; i < u.size(); i++) {
const int M = 20; // output sample points per interval
for (int j = 0; j < M; j++) {
double x = (j + .5) / M;
std::cout << i + x << " " << p[i](x - .5) << " " << u0((i + x) * h - t) << std::endl;
}
std::cout << std::endl;
}
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment