Skip to content
Snippets Groups Projects
Commit 0d313ce8 authored by jstark's avatar jstark
Browse files

Completing 1D unit-test for ENO/WENO scheme.

parent 4a973fdf
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/FiniteDifference/tests/Eno_Weno_unit_test.cpp
+ 147
− 21
View file @ 0d313ce8
......@@ -18,48 +18,144 @@
BOOST_AUTO_TEST_SUITE(EnoWenoTestSuite)
const size_t f_gaussian = 0;
const size_t df_gaussian = 1;
const size_t ddf_gaussian = 2;
const size_t Error1 = 3;
const size_t Error2 = 4;
const size_t ENO_plus = 2;
const size_t ENO_minus = 3;
const size_t Error_plus = 4;
const size_t Error_minus = 5;
BOOST_AUTO_TEST_CASE(Eno_1D_1stDerivative_test)
{
for (size_t N = 64; N <= 64; N *= 2)
for (size_t N = 32; N <= 256; N *= 2)
{
const size_t grid_dim = 1;
const size_t sz[grid_dim] = {N};
const double box_lower = -1.0;
const double box_upper = 1.0;
Box<grid_dim, double> box({box_lower}, {box_upper});
Ghost<grid_dim, long int> ghost(0);
typedef aggregate<double, Point<grid_dim, double>, Point<grid_dim, double>, double, double> props;
Ghost<grid_dim, long int> ghost(3);
typedef aggregate<double, double, double, double, double, double> props;
typedef grid_dist_id<grid_dim, double, props> grid_in_type;
grid_in_type g_dist(sz, box, ghost);
g_dist.setPropNames({"f_gaussian", "df_gaussian", "ddf_gaussian", "Error gradient", "Error Laplace"});
g_dist.setPropNames({"f_gaussian", "df_gaussian", "ENO_plus", "ENO_minus", "Error_plus", "Error_minus"});
double mu = 0.5 * (box_upper - abs(box_lower));
double sigma = 0.1 * (box_upper - box_lower);
auto gdom = g_dist.getDomainGhostIterator();
while (gdom.isNext())
{
auto key = gdom.get();
Point<grid_dim, double> p = g_dist.getPos(key);
// Initialize grid and ghost with gaussian function
g_dist.getProp<f_gaussian>(key) = gaussian(p, mu, sigma);
++gdom;
}
auto dom = g_dist.getDomainIterator();
while (dom.isNext())
{
auto key = dom.get();
Point<grid_dim, double> p = g_dist.getPos(key);
g_dist.getProp<f_gaussian>(key) = gaussian(p, mu, sigma);
for (int d = 0; d < grid_dim; d++)
{
g_dist.getProp<df_gaussian>(key)[d] = gaussian_dx(p.get(d), mu, sigma);
// g_dist.getProp<df_gaussian>(key)[d] = 1;
}
// Analytical solution
g_dist.getProp<df_gaussian>(key) = gaussian(p, mu, sigma, 1);
// ENO_plus
g_dist.getProp<ENO_plus>(key) = ENO_3_Plus<f_gaussian>(g_dist, key, 0);
// ENO_minus
g_dist.getProp<ENO_minus>(key) = ENO_3_Minus<f_gaussian>(g_dist, key, 0);
++dom;
}
// Get the error between analytical and numerical solution
get_relative_error<df_gaussian, ENO_plus, Error_plus>(g_dist);
get_relative_error<df_gaussian, ENO_minus, Error_minus>(g_dist);
{
L_norms lNorms;
lNorms = get_l_norms_grid<Error_plus>(g_dist);
write_lnorms_to_file(N, lNorms, "l_norms_ENO_plus", "./");
}
{
L_norms lNorms;
lNorms = get_l_norms_grid<Error_minus>(g_dist);
write_lnorms_to_file(N, lNorms, "l_norms_ENO_minus", "./");
}
g_dist.write("grid_gaussian_1D_N" + std::to_string(N), FORMAT_BINARY);
g_dist.write("grid_gaussian_ENO_1D_N" + std::to_string(N), FORMAT_BINARY);
}
}
const size_t WENO_plus = 2;
const size_t WENO_minus = 3;
BOOST_AUTO_TEST_CASE(Weno_1D_1stDerivative_test)
{
for (size_t N = 32; N <= 256; N *= 2)
{
const size_t grid_dim = 1;
const size_t sz[grid_dim] = {N};
const double box_lower = -1.0;
const double box_upper = 1.0;
Box<grid_dim, double> box({box_lower}, {box_upper});
Ghost<grid_dim, long int> ghost(3);
typedef aggregate<double, double, double, double, double, double> props;
typedef grid_dist_id<grid_dim, double, props> grid_in_type;
grid_in_type g_dist(sz, box, ghost);
g_dist.setPropNames({"f_gaussian", "df_gaussian", "WENO_plus", "WENO_minus", "Error_plus", "Error_minus"});
double mu = 0.5 * (box_upper - abs(box_lower));
double sigma = 0.1 * (box_upper - box_lower);
auto gdom = g_dist.getDomainGhostIterator();
while (gdom.isNext())
{
auto key = gdom.get();
Point<grid_dim, double> p = g_dist.getPos(key);
// Initialize grid and ghost with gaussian function
g_dist.getProp<f_gaussian>(key) = gaussian(p, mu, sigma);
++gdom;
}
auto dom = g_dist.getDomainIterator();
while (dom.isNext())
{
auto key = dom.get();
Point<grid_dim, double> p = g_dist.getPos(key);
// Analytical solution
g_dist.getProp<f_gaussian>(key) = gaussian(p, mu, sigma);
g_dist.getProp<df_gaussian>(key) = gaussian(p, mu, sigma, 1);
// ENO_plus
g_dist.getProp<WENO_plus>(key) = WENO_5_Plus<f_gaussian>(g_dist, key, 0);
// ENO_minus
g_dist.getProp<WENO_minus>(key) = WENO_5_Minus<f_gaussian>(g_dist, key, 0);
++dom;
}
// Get the error between analytical and numerical solution
get_relative_error<df_gaussian, WENO_plus, Error_plus>(g_dist);
get_relative_error<df_gaussian, WENO_minus, Error_minus>(g_dist);
{
L_norms lNorms;
lNorms = get_l_norms_grid<Error_plus>(g_dist);
write_lnorms_to_file(N, lNorms, "l_norms_WENO_plus", "./");
}
{
L_norms lNorms;
lNorms = get_l_norms_grid<Error_minus>(g_dist);
write_lnorms_to_file(N, lNorms, "l_norms_WENO_minus", "./");
}
g_dist.write("grid_gaussian_WENO_1D_N" + std::to_string(N), FORMAT_BINARY);
}
}
BOOST_AUTO_TEST_CASE(Eno_2D_1stDerivative_test)
{
for(size_t N=32; N<=128; N*=2)
......@@ -69,28 +165,58 @@ BOOST_AUTO_TEST_SUITE(EnoWenoTestSuite)
const double box_lower = -1.0;
const double box_upper = 1.0;
Box<grid_dim, double> box({box_lower, box_lower}, {box_upper, box_upper});
Ghost<grid_dim, long int> ghost(0);
typedef aggregate<double, Point<grid_dim, double>, Point<grid_dim, double>, double, double> props;
Ghost<grid_dim, long int> ghost(3);
typedef aggregate<double, Point<grid_dim, double>, Point<grid_dim, double>, Point<grid_dim, double>, double, double> props;
typedef grid_dist_id<grid_dim, double, props> grid_in_type;
grid_in_type g_dist(sz, box, ghost);
g_dist.setPropNames({"f_gaussian", "df_gaussian", "ddf_gaussian", "Error gradient", "Error Laplace"});
g_dist.setPropNames({"f_gaussian", "df_gaussian", "ENO_plus", "ENO_minus", "Error_plus", "Error_minus"});
double mu = 0.5 * (box_upper - abs(box_lower));
double sigma = 0.1 * (box_upper - box_lower);
auto gdom = g_dist.getDomainGhostIterator();
while (gdom.isNext())
{
auto key = gdom.get();
Point<grid_dim, double> p = g_dist.getPos(key);
// Initialize grid and ghost with gaussian function
g_dist.getProp<f_gaussian>(key) = gaussian(p, mu, sigma);
++gdom;
}
auto dom = g_dist.getDomainIterator();
while (dom.isNext())
{
auto key = dom.get();
Point<grid_dim, double> p = g_dist.getPos(key);
g_dist.getProp<f_gaussian>(key) = gaussian(p, mu, sigma);
for (int d = 0; d < grid_dim; d++)
{
g_dist.getProp<df_gaussian>(key)[d] = gaussian_dx(p.get(d), mu, sigma);
g_dist.getProp<df_gaussian>(key)[d] = gaussian(p.get(d), mu, sigma, 1);
g_dist.getProp<ENO_plus>(key)[d] = ENO_3_Plus<f_gaussian>(g_dist, key, d);
g_dist.getProp<ENO_minus>(key)[d] = ENO_3_Minus<f_gaussian>(g_dist, key, d);
}
++dom;
}
// Get the error between analytical and numerical solution
get_relative_error<df_gaussian, ENO_plus, Error_plus>(g_dist);
get_relative_error<df_gaussian, ENO_minus, Error_minus>(g_dist);
{
L_norms lNorms;
lNorms = get_l_norms_grid<Error_plus>(g_dist);
write_lnorms_to_file(N, lNorms, "l_norms_2D_ENO_plus", "./");
}
{
L_norms lNorms;
lNorms = get_l_norms_grid<Error_minus>(g_dist);
write_lnorms_to_file(N, lNorms, "l_norms_2D_ENO_minus", "./");
}
g_dist.write("grid_gaussian_2D_N" + std::to_string(N), FORMAT_BINARY);
}
}
......
src/FiniteDifference/tests/Gaussian.hpp
+ 35
− 44
View file @ 0d313ce8
......@@ -7,71 +7,62 @@
#include "cmath"
double hermite_polynomial(double x, double sigma, int order)
{
double h;
switch(order)
{
case 0:
h = 1;
break;
case 1:
h = -x / (sigma * sigma);
break;
case 2:
h = (x*x - sigma*sigma) / (sigma*sigma*sigma*sigma);
break;
default:
h = 0;
std::cout << "Only gaussian derivatives of order 0, 1, and 2 implemented. Aborting..." << std::endl;
abort();
break;
}
return h;
}
template <typename point_type>
double gaussian(const point_type & x, const double mu, const double sigma)
double gaussian(const point_type & x, const double mu, const double sigma, const int order=0)
{
const double pi = 3.14159265358979323846;
const double sqrt2pi = sqrt(2*pi);
const double sigmapow2 = sigma * sigma;
double sum = 0;
double normalization_factor = 1;
double h = 1;
for(int d=0; d<point_type::dims; d++)
{
sum += (x.get(d) - mu) * (x.get(d) - mu) / (sigmapow2);;
sum += (x.get(d) - mu) * (x.get(d) - mu) / (sigma*sigma);;
normalization_factor *= 1 / (sqrt2pi * sigma);
h *= hermite_polynomial(x.get(d), sigma, order);
}
return normalization_factor * exp(-0.5 * sum);
return h * normalization_factor * exp(-0.5 * sum);
}
// 1D
double gaussian(const double x, const double mu, const double sigma)
double gaussian(const double x, const double mu, const double sigma, const int order=0)
{
const double pi = 3.14159265358979323846;
const double sqrt2pi = sqrt(2*pi);
const double sigmapow2 = sigma * sigma;
double sum = 0;
double normalization_factor = 1;
double sum = (x - mu) * (x - mu) / (sigma*sigma);;
double normalization_factor = 1 / (sqrt2pi * sigma);
double h = hermite_polynomial(x, sigma, order);
sum += (x - mu) * (x - mu) / (sigmapow2);;
normalization_factor *= 1 / (sqrt2pi * sigma);
return normalization_factor * exp(-0.5 * sum);
return h * normalization_factor * exp(-0.5 * sum);
}
template <typename point_type>
double gaussian_dx(const point_type & x, const double mu, const double sigma)
{
const double pi = 3.14159265358979323846;
const double sqrt2pi = sqrt(2*pi);
const double sigmapow2 = sigma * sigma;
double sum = (x - mu) * (x - mu) / (sigmapow2);
double normalization_factor = 1 /( sqrt2pi * sigma);
double prefactor = - x / sigmapow2;
return prefactor * normalization_factor * exp(-0.5 * sum);
}
template <typename point_type>
double gaussian_ddx(const point_type & x, const double mu, const double sigma)
{
const double pi = 3.14159265358979323846;
const double sqrt2pi = sqrt(2*pi);
const double sigmapow2 = sigma * sigma;
const double sigmapow4 = sigmapow2 * sigmapow2;
double sum = (x - mu) * (x - mu) / (sigmapow2);
double normalization_factor = 1 /( sqrt2pi * sigma);
double prefactor = - (sigmapow2 - x*x) / sigmapow4;
return prefactor * normalization_factor * exp(-0.5 * sum);
}
#endif //OPENFPM_NUMERICS_GAUSSIAN_HPP
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment