Skip to content
Snippets Groups Projects
Commit bbe9e090 authored by Abhinav Singh's avatar Abhinav Singh
Browse files

Merge remote-tracking branch 'origin/develop' into develop

parents f7094e9e 8911e4b7
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/FiniteDifference/Upwind_gradient.hpp
+ 44
− 19
View file @ bbe9e090
......@@ -27,12 +27,13 @@
#include "Grid/grid_dist_id.hpp"
#include "FD_simple.hpp"
#include "Eno_Weno.hpp"
#include "level_set/redistancing_Sussman/HelpFunctions.hpp"
/**@brief Upwinding: For a specific dimension, from the forward and backward gradient find the upwind side.
*
* @param dplus: Gradient approximated using RHS neighbors.
* @param dminus: Gradient approximated using LHS neighbors.
* @param sign: Sign of the velocity with which the wave front is moving.
* @param sign: Sign of current component of the wave front velocity.
* @return Scalar upwind gradient approximation in the dimension given.
*/
template <typename field_type>
......@@ -50,14 +51,38 @@ static field_type upwinding(field_type dplus, field_type dminus, int sign)
return grad_upwind;
}
/**@brief Returns the sign of a scalar velocity field
*
* @tparam velocity_type Template type of velocity (scalar).
* @param v Velocity of type velocity_type.
* @param d Dimension. This is a dummy to simplify call in module #FD_upwind.
* @return Sign of velocity #v.
*/
template <typename velocity_type>
inline int get_sign_velocity(velocity_type v, size_t d)
{
return sgn(v);
}
/**@brief Returns the sign of one component in a vector velocity field
*
* @tparam velocity_type Template type of velocity (vector).
* @param v Velocity of type velocity_type.
* @param d Component of velocity vector for which sign should be returned.
* @return Sign of component d from velocity #v.
*/
template <typename velocity_type>
inline int get_sign_velocity(velocity_type v[], size_t d)
{
return sgn(v[d]);
}
/**@brief Get the upwind finite difference of a scalar property on the current grid node.
*
* @details Order of accuracy can be 1, 3 or 5.
*
* @tparam Field Size_t index of property for which the gradient should be computed.
* @tparam Sign Size_t index of property that contains the initial sign of that property for which the upwind FD
* should be computed.
* @tparam Velocity Size_t index of property that contains the velocity field. Can be scalar or vector field.
* @tparam gridtype Type of input grid.
* @tparam keytype Type of key variable.
* @param grid Grid, on which the gradient should be computed.
......@@ -66,14 +91,16 @@ static field_type upwinding(field_type dplus, field_type dminus, int sign)
* @param order Order of accuracy of the difference scheme. Can be 1, 3 or 5.
* @return Upwind finite difference in one dimension of the property under index Field on the current node with index key.
*/
template <size_t Field, size_t Sign, typename gridtype, typename keytype>
template <size_t Field, size_t Velocity, typename gridtype, typename keytype>
auto FD_upwind(gridtype &grid, keytype &key, size_t d, size_t order)
{
// std::cout << boost::typeindex::type_id_with_cvr<std::decay_t<decltype(first_node)>>() << std::endl;
typedef typename std::decay_t<decltype(grid.template get<Field>(key))> field_type;
field_type dplus = 0.0, dminus = 0.0;
int sign_phi0 = grid.template get<Sign> (key);
// Get the sign of the velocity for the current dimension
int sign_velocity = get_sign_velocity(grid.template get<Field>(key), d);
switch(order)
{
......@@ -98,7 +125,7 @@ auto FD_upwind(gridtype &grid, keytype &key, size_t d, size_t order)
break;
}
return upwinding(dplus, dminus, sign_phi0);
return upwinding(dplus, dminus, sign_velocity);
}
......@@ -109,8 +136,7 @@ auto FD_upwind(gridtype &grid, keytype &key, size_t d, size_t order)
* respectively, depending on the side of the border.
*
* @tparam Field Size_t index of property for which the gradient should be computed.
* @tparam Sign Size_t index of property that contains the initial sign of that property for which the upwind FD
* should be computed.
* @tparam Velocity Size_t index of property that contains the velocity field. Can be scalar or vector field.
* @tparam Gradient Size_t index of property where the gradient result should be stored.
* @tparam gridtype Type of input grid.
* @param grid Grid, on which the gradient should be computed.
......@@ -119,11 +145,11 @@ auto FD_upwind(gridtype &grid, keytype &key, size_t d, size_t order)
* @param order Size_t variable, order of accuracy the upwind FD scheme should have. Can be 1, 3 or 5.
*/
// Use upwinding for inner grid points and one sided backward / forward stencil at border (if one_sided_BC=true)
template <size_t Field, size_t Sign, size_t Gradient, typename gridtype>
template <size_t Field, size_t Velocity, size_t Gradient, typename gridtype>
void upwind_gradient(gridtype & grid, const bool one_sided_BC, size_t order)
{
grid.template ghost_get<Field>(KEEP_PROPERTIES);
grid.template ghost_get<Sign>(KEEP_PROPERTIES);
grid.template ghost_get<Velocity>(KEEP_PROPERTIES);
auto dom = grid.getDomainIterator();
......@@ -140,13 +166,13 @@ void upwind_gradient(gridtype & grid, const bool one_sided_BC, size_t order)
if (key_g.get(d) > 2 && key_g.get(d) < grid.size(d) - 3) // if point lays with min. 3 nodes distance to
// boundary
{
grid.template get<Gradient> (key) [d] = FD_upwind<Field, Sign>(grid, key, d, order);
grid.template get<Gradient> (key) [d] = FD_upwind<Field, Velocity>(grid, key, d, order);
}
// Grid nodes in stencil-wide vicinity to boundary
else if (key_g.get(d) > 0 && key_g.get(d) < grid.size(d) - 1) // if point lays not on the grid boundary
{
grid.template get<Gradient> (key) [d] = FD_upwind<Field, Sign>(grid, key, d, 1);
grid.template get<Gradient> (key) [d] = FD_upwind<Field, Velocity>(grid, key, d, 1);
}
else if (key_g.get(d) == 0) // if point lays at left boundary, use right sided kernel
......@@ -170,7 +196,7 @@ void upwind_gradient(gridtype & grid, const bool one_sided_BC, size_t order)
auto key = dom.get();
for(size_t d = 0; d < gridtype::dims; d++)
{
grid.template get<Gradient> (key) [d] = FD_upwind<Field, Sign>(grid, key, d, order);
grid.template get<Gradient> (key) [d] = FD_upwind<Field, Velocity>(grid, key, d, order);
}
++dom;
}
......@@ -207,17 +233,16 @@ static bool ghost_width_is_sufficient(gridtype & grid, size_t required_width)
/**@brief Calls #upwind_gradient. Computes upwind gradient of desired order {1, 3, 5} for the whole n-dim grid.
*
* @tparam Field Size_t index of property for which the gradient should be computed.
* @tparam Sign Size_t index of property that contains the initial sign of that property for which the upwind FD
* should be computed / sign of the velocity.
* @tparam Velocity Size_t index of property that contains the velocity field. Can be scalar or vector field.
* @tparam Gradient Size_t index of property where the upwind gradient result should be stored.
* @tparam gridtype Type of input grid.
* @param grid Grid, on which the gradient should be computed.
*/
template <size_t Field_in, size_t Sign, size_t Gradient_out, typename gridtype>
template <size_t Field_in, size_t Velocity, size_t Gradient_out, typename gridtype>
void get_upwind_gradient(gridtype & grid, const size_t order=5, const bool one_sided_BC=true)
{
grid.template ghost_get<Field_in>();
grid.template ghost_get<Sign>();
grid.template ghost_get<Velocity>();
if (!one_sided_BC)
{
......@@ -241,13 +266,13 @@ void get_upwind_gradient(gridtype & grid, const size_t order=5, const bool one_s
case 1:
case 3:
case 5:
upwind_gradient<Field_in, Sign, Gradient_out>(grid, one_sided_BC, order);
upwind_gradient<Field_in, Velocity, Gradient_out>(grid, one_sided_BC, order);
break;
default:
auto &v_cl = create_vcluster();
if (v_cl.rank() == 0) std::cout << "Order of accuracy chosen not valid. Using default order 1." <<
std::endl;
upwind_gradient<Field_in, Sign, Gradient_out>(grid, one_sided_BC, 1);
upwind_gradient<Field_in, Velocity, Gradient_out>(grid, one_sided_BC, 1);
break;
}
......
src/FiniteDifference/tests/Upwind_gradient_unit_test.cpp
+ 123
− 33
View file @ bbe9e090
......@@ -13,11 +13,11 @@
BOOST_AUTO_TEST_SUITE(UpwindGradientTestSuite)
const double EPSILON = 1e-14;
const size_t f_gaussian = 0;
const size_t Sign = 1;
const size_t df_gaussian = 2;
const size_t df_upwind = 3;
const size_t Error = 4;
const size_t F_GAUSSIAN = 0;
const size_t VELOCITY = 1;
const size_t dF_GAUSSIAN = 2;
const size_t dF_UPWIND = 3;
const size_t ERROR = 4;
BOOST_AUTO_TEST_CASE(Upwind_gradient_order1_1D_test)
{
......@@ -40,7 +40,7 @@ BOOST_AUTO_TEST_SUITE(UpwindGradientTestSuite)
const size_t sz[grid_dim] = {N};
grid_in_type g_dist(sz, box, ghost);
g_dist.setPropNames({"f_gaussian", "Sign", "df_gaussian", "df_upwind", "Error"});
g_dist.setPropNames({"F_GAUSSIAN", "VELOCITY", "dF_GAUSSIAN", "dF_UPWIND", "ERROR"});
auto gdom = g_dist.getDomainGhostIterator();
while (gdom.isNext())
......@@ -48,9 +48,9 @@ BOOST_AUTO_TEST_SUITE(UpwindGradientTestSuite)
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);
g_dist.getProp<F_GAUSSIAN>(key) = gaussian(p, mu, sigma);
// Initialize the sign of the gaussian function
g_dist.getProp<Sign>(key) = sgn(g_dist.getProp<f_gaussian>(key));
g_dist.getProp<VELOCITY>(key) = sgn(g_dist.getProp<F_GAUSSIAN>(key));
++gdom;
}
......@@ -63,19 +63,19 @@ BOOST_AUTO_TEST_SUITE(UpwindGradientTestSuite)
for (int d = 0; d < grid_dim; d++)
{
// Analytical gradient
g_dist.getProp<df_gaussian>(key)[d] = hermite_polynomial(p.get(d), sigma, 1) * g_dist.getProp<f_gaussian>(key);
g_dist.getProp<dF_GAUSSIAN>(key)[d] = hermite_polynomial(p.get(d), sigma, 1) * g_dist.getProp<F_GAUSSIAN>(key);
}
++dom;
}
// Upwind gradient with specific convergence order, not becoming one-sided at the boundary
get_upwind_gradient<f_gaussian, Sign, df_upwind>(g_dist, convergence_order, false);
get_upwind_gradient<F_GAUSSIAN, VELOCITY, dF_UPWIND>(g_dist, convergence_order, false);
// Get the error between analytical and numerical solution
get_relative_error<df_upwind, df_gaussian, Error>(g_dist);
get_relative_error<dF_UPWIND, dF_GAUSSIAN, ERROR>(g_dist);
LNorms<double> lNorms;
lNorms.get_l_norms_grid<Error>(g_dist);
lNorms.get_l_norms_grid<ERROR>(g_dist);
if (N==32) BOOST_CHECK_MESSAGE(lNorms.l2 <= 0.38954184748195 + EPSILON, "Checking L2-norm upwind gradient "
"order " + std::to_string
......@@ -105,7 +105,7 @@ BOOST_AUTO_TEST_SUITE(UpwindGradientTestSuite)
const size_t sz[grid_dim] = {N};
grid_in_type g_dist(sz, box, ghost);
g_dist.setPropNames({"f_gaussian", "Sign", "df_gaussian", "df_upwind", "Error"});
g_dist.setPropNames({"F_GAUSSIAN", "VELOCITY", "dF_GAUSSIAN", "dF_UPWIND", "ERROR"});
auto gdom = g_dist.getDomainGhostIterator();
while (gdom.isNext())
......@@ -113,9 +113,9 @@ BOOST_AUTO_TEST_SUITE(UpwindGradientTestSuite)
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);
g_dist.getProp<F_GAUSSIAN>(key) = gaussian(p, mu, sigma);
// Initialize the sign of the gaussian function
g_dist.getProp<Sign>(key) = sgn(g_dist.getProp<f_gaussian>(key));
g_dist.getProp<VELOCITY>(key) = sgn(g_dist.getProp<F_GAUSSIAN>(key));
++gdom;
}
......@@ -128,19 +128,19 @@ BOOST_AUTO_TEST_SUITE(UpwindGradientTestSuite)
for (int d = 0; d < grid_dim; d++)
{
// Analytical gradient
g_dist.getProp<df_gaussian>(key)[d] = hermite_polynomial(p.get(d), sigma, 1) * g_dist.getProp<f_gaussian>(key);
g_dist.getProp<dF_GAUSSIAN>(key)[d] = hermite_polynomial(p.get(d), sigma, 1) * g_dist.getProp<F_GAUSSIAN>(key);
}
++dom;
}
// Upwind gradient with specific convergence order, not becoming one-sided at the boundary
get_upwind_gradient<f_gaussian, Sign, df_upwind>(g_dist, convergence_order, false);
get_upwind_gradient<F_GAUSSIAN, VELOCITY, dF_UPWIND>(g_dist, convergence_order, false);
// Get the error between analytical and numerical solution
get_relative_error<df_upwind, df_gaussian, Error>(g_dist);
get_relative_error<dF_UPWIND, dF_GAUSSIAN, ERROR>(g_dist);
LNorms<double> lNorms;
lNorms.get_l_norms_grid<Error>(g_dist);
lNorms.get_l_norms_grid<ERROR>(g_dist);
if (N==32) BOOST_CHECK_MESSAGE(lNorms.l2 <= 0.08667855716144 + EPSILON, "Checking L2-norm upwind gradient "
"order "
......@@ -168,7 +168,7 @@ BOOST_AUTO_TEST_SUITE(UpwindGradientTestSuite)
const size_t sz[grid_dim] = {N};
grid_in_type g_dist(sz, box, ghost);
g_dist.setPropNames({"f_gaussian", "Sign", "df_gaussian", "df_upwind", "Error"});
g_dist.setPropNames({"F_GAUSSIAN", "VELOCITY", "dF_GAUSSIAN", "dF_UPWIND", "ERROR"});
auto gdom = g_dist.getDomainGhostIterator();
while (gdom.isNext())
......@@ -176,9 +176,9 @@ BOOST_AUTO_TEST_SUITE(UpwindGradientTestSuite)
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);
g_dist.getProp<F_GAUSSIAN>(key) = gaussian(p, mu, sigma);
// Initialize the sign of the gaussian function
g_dist.getProp<Sign>(key) = sgn(g_dist.getProp<f_gaussian>(key));
g_dist.getProp<VELOCITY>(key) = sgn(g_dist.getProp<F_GAUSSIAN>(key));
++gdom;
}
......@@ -191,19 +191,19 @@ BOOST_AUTO_TEST_SUITE(UpwindGradientTestSuite)
for (int d = 0; d < grid_dim; d++)
{
// Analytical gradient
g_dist.getProp<df_gaussian>(key)[d] = hermite_polynomial(p.get(d), sigma, 1) * g_dist.getProp<f_gaussian>(key);
g_dist.getProp<dF_GAUSSIAN>(key)[d] = hermite_polynomial(p.get(d), sigma, 1) * g_dist.getProp<F_GAUSSIAN>(key);
}
++dom;
}
// Upwind gradient with specific convergence order, not becoming one-sided at the boundary
get_upwind_gradient<f_gaussian, Sign, df_upwind>(g_dist, convergence_order, false);
get_upwind_gradient<F_GAUSSIAN, VELOCITY, dF_UPWIND>(g_dist, convergence_order, false);
// Get the error between analytical and numerical solution
get_relative_error<df_upwind, df_gaussian, Error>(g_dist);
get_relative_error<dF_UPWIND, dF_GAUSSIAN, ERROR>(g_dist);
LNorms<double> lNorms;
lNorms.get_l_norms_grid<Error>(g_dist);
lNorms.get_l_norms_grid<ERROR>(g_dist);
if (N==32) BOOST_CHECK_MESSAGE(lNorms.l2 <= 0.03215172234342 + EPSILON, "Checking L2-norm upwind gradient "
"order "
......@@ -229,7 +229,7 @@ BOOST_AUTO_TEST_SUITE(UpwindGradientTestSuite)
{
const size_t sz[grid_dim] = {N, N, N};
grid_in_type g_dist(sz, box, ghost);
g_dist.setPropNames({"f_gaussian", "Sign", "df_gaussian", "df_upwind", "Error"});
g_dist.setPropNames({"F_GAUSSIAN", "VELOCITY", "dF_GAUSSIAN", "dF_UPWIND", "ERROR"});
auto gdom = g_dist.getDomainGhostIterator();
while (gdom.isNext())
......@@ -237,9 +237,9 @@ BOOST_AUTO_TEST_SUITE(UpwindGradientTestSuite)
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);
g_dist.getProp<F_GAUSSIAN>(key) = gaussian(p, mu, sigma);
// Initialize the sign of the gaussian function
g_dist.getProp<Sign>(key) = sgn(g_dist.getProp<f_gaussian>(key));
g_dist.getProp<VELOCITY>(key) = sgn(g_dist.getProp<F_GAUSSIAN>(key));
++gdom;
}
......@@ -252,7 +252,7 @@ BOOST_AUTO_TEST_SUITE(UpwindGradientTestSuite)
for (int d = 0; d < grid_dim; d++)
{
// Analytical gradient
g_dist.getProp<df_gaussian>(key)[d] = hermite_polynomial(p.get(d), sigma, 1) * g_dist.getProp<f_gaussian>(key);
g_dist.getProp<dF_GAUSSIAN>(key)[d] = hermite_polynomial(p.get(d), sigma, 1) * g_dist.getProp<F_GAUSSIAN>(key);
}
++dom;
}
......@@ -260,12 +260,12 @@ BOOST_AUTO_TEST_SUITE(UpwindGradientTestSuite)
for (int convergence_order = 1; convergence_order <=5; convergence_order +=2 )
{
// Upwind gradient with specific convergence order, not becoming one-sided at the boundary
get_upwind_gradient<f_gaussian, Sign, df_upwind>(g_dist, convergence_order, false);
get_upwind_gradient<F_GAUSSIAN, VELOCITY, dF_UPWIND>(g_dist, convergence_order, false);
// Get the error between analytical and numerical solution
get_relative_error<df_upwind, df_gaussian, Error>(g_dist);
get_relative_error<dF_UPWIND, dF_GAUSSIAN, ERROR>(g_dist);
LNorms<double> lNorms;
lNorms.get_l_norms_grid<Error>(g_dist);
lNorms.get_l_norms_grid<ERROR>(g_dist);
if(N==32)
{
......@@ -290,4 +290,94 @@ BOOST_AUTO_TEST_SUITE(UpwindGradientTestSuite)
}
}
}
BOOST_AUTO_TEST_CASE(Upwind_gradient_3D_vector_velocity_test)
{
{
const size_t grid_dim = 3;
const double box_lower = -1.0;
const double box_upper = 1.0;
Box<grid_dim, double> box({box_lower, box_lower, box_lower}, {box_upper, box_upper, box_upper});
Ghost<grid_dim, long int> ghost(3);
typedef aggregate<double, double[grid_dim], Point<grid_dim, double>, Point<grid_dim, double>, double> props;
typedef grid_dist_id<grid_dim, double, props> grid_in_type;
double mu = 0.5 * (box_upper - abs(box_lower));
double sigma = 0.3 * (box_upper - box_lower);
size_t N = 32;
// for(size_t N = 32; N <=256; N *=2)
{
const size_t sz[grid_dim] = {N, N, N};
grid_in_type g_dist(sz, box, ghost);
g_dist.setPropNames({"F_GAUSSIAN", "Velocity", "dF_GAUSSIAN", "dF_UPWIND", "ERROR"});
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);
for (int d = 0; d < grid_dim; d++)
{
// Analytical gradient
g_dist.getProp<dF_GAUSSIAN>(key)[d] = hermite_polynomial(p.get(d), sigma, 1) * g_dist.getProp<F_GAUSSIAN>(key);
// Initialize the velocity vector field
g_dist.getProp<VELOCITY>(key)[d] = g_dist.getProp<dF_GAUSSIAN>(key)[d];
}
++dom;
}
for (int convergence_order = 1; convergence_order <=5; convergence_order +=2 )
{
// Upwind gradient with specific convergence order, not becoming one-sided at the boundary
get_upwind_gradient<F_GAUSSIAN, VELOCITY, dF_UPWIND>(g_dist, convergence_order, false);
// Get the error between analytical and numerical solution
get_relative_error<dF_UPWIND, dF_GAUSSIAN, ERROR>(g_dist);
LNorms<double> lNorms;
lNorms.get_l_norms_grid<ERROR>(g_dist);
if(N==32)
{
switch(convergence_order)
{
case 1:
BOOST_CHECK_MESSAGE(lNorms.l2 <= 0.38954184748194 + EPSILON, "Checking L2-norm upwind"
" gradient with vector "
"velocity order " +
std::to_string(convergence_order));
break;
case 3:
BOOST_CHECK_MESSAGE(lNorms.l2 <= 0.08667855716144 + EPSILON, "Checking L2-norm upwind"
" gradient with "
"vector velocity order " +
std::to_string(convergence_order));
break;
case 5:
BOOST_CHECK_MESSAGE(lNorms.l2 <= 0.02285996528578 + EPSILON, "Checking L2-norm upwind"
" gradient with "
"vector velocity order " +
std::to_string(convergence_order));
break;
default:
std::cout << "Checking only implemented for convergence order 1, 3 and 5." << std::endl;
}
}
}
}
}
}
BOOST_AUTO_TEST_SUITE_END()
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