Skip to content

Commit

Permalink
Add test
Browse files Browse the repository at this point in the history
  • Loading branch information
@gassmoeller
gassmoeller committed Apr 26, 2018
1 parent 402da57 commit 5a8e038
Show file tree
Hide file tree
Showing 3 changed files with 205 additions and 0 deletions.
95 changes: 95 additions & 0 deletions tests/ellipsoidal_chunk_geometry.cc
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,95 @@
#include <aspect/geometry_model/ellipsoidal_chunk.h>
#include <aspect/geometry_model/initial_topography_model/zero_topography.h>

#include <deal.II/base/exceptions.h>

#include <iostream>

using namespace aspect;
using namespace dealii;

bool test_point(const GeometryModel::EllipsoidalChunk<3>::EllipsoidalChunkGeometry ellipsoidal_manifold,
const Point<3> &test_point)
{
const Point<3> converted_point = ellipsoidal_manifold.pull_back(test_point);
const Point<3> twice_converted_point = ellipsoidal_manifold.push_forward(converted_point);

const Tensor<1,3> difference = twice_converted_point - test_point;
const double error = difference.norm() / test_point.norm();

std::cout << "Point: " << test_point
<< " becomes: " << converted_point
<< ". Back transformed: " << twice_converted_point
<< ". Error: " << error << std::endl;

return true;
}

/*
* Launch the following function when this plugin is created.
*/
int f()
{

const unsigned int dim=3;

InitialTopographyModel::ZeroTopography<dim> topography;
GeometryModel::EllipsoidalChunk<dim>::EllipsoidalChunkGeometry ellipsoidal_manifold;
ellipsoidal_manifold.initialize(&topography);

std::vector<Point<2> > corners(2,Point<2>(-15.0,-15.0));
corners[1] *= -1.0;

std::cout << "Simple sphere test" << std::endl;
ellipsoidal_manifold.set_manifold_parameters(6371000.0,
0.0,
6371000.0,
2890000.0,
corners);

std::vector<Point<3> > test_points;
test_points.push_back(Point<3> (6371000.0,0,0));
test_points.push_back(Point<3> (6171000.0,0,0));
test_points.push_back(Point<3> (3000000.0,3000000.0,0));
test_points.push_back(Point<3> (3000000.0,3000000.0,3000000.0));
test_points.push_back(Point<3> (3000000.0,-3000000.0,-3000000.0));
test_points.push_back(Point<3> (25000.0,25000.0,3481000.0));
test_points.push_back(Point<3> (25000.0,25000.0,1000000.0));
test_points.push_back(Point<3> (25000.0,25000.0,500000.0));
test_points.push_back(Point<3> (25000.0,25000.0,50000.0));

for (unsigned int i=0; i<test_points.size(); ++i)
test_point(ellipsoidal_manifold,test_points[i]);

std::cout << "WGS84 test" << std::endl;
const double semi_major_axis_a = 6378137.0;
const double eccentricity = 8.1819190842622e-2;
const double semi_minor = std::sqrt((1 - pow(eccentricity,2)) * pow(semi_major_axis_a,2));

ellipsoidal_manifold.set_manifold_parameters(semi_major_axis_a,
eccentricity,
semi_minor,
2890000.0,
corners);

for (unsigned int i=0; i<test_points.size(); ++i)
test_point(ellipsoidal_manifold,test_points[i]);


std::cout << "Test for points outside of the provided depth range" << std::endl;
ellipsoidal_manifold.set_manifold_parameters(semi_major_axis_a,
eccentricity,
semi_minor,
500000.0,
corners);

for (unsigned int i=0; i<test_points.size(); ++i)
test_point(ellipsoidal_manifold,test_points[i]);

exit (0);
return 42;
}


// run this function by initializing a global variable by it
int i = f();
76 changes: 76 additions & 0 deletions tests/ellipsoidal_chunk_geometry.prm
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,76 @@
# This test checks the accuracy of the coordinate transformation
# for the ellipsoidal chunk geometry model. Because the used
# manifold includes an approximation for the conversion between
# spherical and ellipsoidal coordinates the conversion is inaccurate
# very close to the center of the domain. This test documents this
# inaccuracy and shows that it is negligible for points at the
# core-mantle boundary (5e-11 relative to center distance) and even for points
# at the inner-core-boundary (1e-7).

set Dimension = 3
set Use years in output instead of seconds = true
set End time = 0

subsection Geometry model
set Model name = ellipsoidal chunk
subsection Ellipsoidal chunk
set NE corner = 7.5:30
set SW corner = -7.5:-10
set Depth = 3000000
set Eccentricity = 0.5
set East-West subdivisions = 1
set North-South subdivisions = 2
set Depth subdivisions = 3
end
end

subsection Boundary velocity model
set Zero velocity boundary indicators = inner
set Tangential velocity boundary indicators = outer,north,south,east,west
end

subsection Boundary temperature model
set Fixed temperature boundary indicators = inner, outer,north,south,east,west
set List of model names = constant
subsection Constant
set Boundary indicator to temperature mappings = west: 0, east: 1, south: 2, north: 3, inner: 4, outer:5
end
end


subsection Material model
set Model name = simple

subsection Simple model
set Thermal expansion coefficient = 4e-5
set Viscosity = 1e22
end
end

subsection Initial temperature model
set Model name = function
subsection Function
set Function expression = 1.473e3
end
end

subsection Gravity model
set Model name = radial constant
end

subsection Mesh refinement
set Initial global refinement = 1
set Initial adaptive refinement = 1
set Strategy = temperature
set Time steps between mesh refinement = 15
end

subsection Postprocess
set List of postprocessors = visualization, velocity statistics, temperature statistics, basic statistics

subsection Visualization
set List of output variables = density
set Output format = gnuplot

end
end
34 changes: 34 additions & 0 deletions tests/ellipsoidal_chunk_geometry/screen-output
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,34 @@
-----------------------------------------------------------------------------
-----------------------------------------------------------------------------

Loading shared library <./libellipsoidal_chunk_geometry.so>
Simple sphere test
Point: 6.371e+06 0 0 becomes: 0 0 0. Back transformed: 6.371e+06 0 0. Error: 0
Point: 6.171e+06 0 0 becomes: 0 0 -200000. Back transformed: 6.171e+06 0 0. Error: 0
Point: 3e+06 3e+06 0 becomes: 0.785398 0 -2.12836e+06. Back transformed: 3e+06 3e+06 0. Error: 1.09757e-16
Point: 3e+06 3e+06 3e+06 becomes: 0.785398 0.61548 -1.17485e+06. Back transformed: 3e+06 3e+06 3e+06. Error: 0
Point: 3e+06 -3e+06 -3e+06 becomes: -0.785398 -0.61548 -1.17485e+06. Back transformed: 3e+06 -3e+06 -3e+06. Error: 0
Point: 25000 25000 3.481e+06 becomes: 0.785398 1.56064 -2.88982e+06. Back transformed: 25000 25000 3.481e+06. Error: 2.14025e-15
Point: 25000 25000 1e+06 becomes: 0.785398 1.53546 -5.37038e+06. Back transformed: 25000 25000 1e+06. Error: 6.98103e-16
Point: 25000 25000 500000 becomes: 0.785398 1.5002 -5.86975e+06. Back transformed: 25000 25000 500000. Error: 9.30731e-16
Point: 25000 25000 50000 becomes: 0.785398 0.955317 -6.30976e+06. Back transformed: 25000 25000 50000. Error: 7.44593e-15
WGS84 test
Point: 6.371e+06 0 0 becomes: 0 0 -7137. Back transformed: 6.371e+06 0 0. Error: 0
Point: 6.171e+06 0 0 becomes: 0 0 -207137. Back transformed: 6.171e+06 0 0. Error: 0
Point: 3e+06 3e+06 0 becomes: 0.785398 0 -2.1355e+06. Back transformed: 3e+06 3e+06 0. Error: 1.09757e-16
Point: 3e+06 3e+06 3e+06 becomes: 0.785398 0.619368 -1.17483e+06. Back transformed: 3e+06 3e+06 3e+06. Error: 3.69023e-09
Point: 3e+06 -3e+06 -3e+06 becomes: -0.785398 -0.619368 -1.17483e+06. Back transformed: 3e+06 -3e+06 -3e+06. Error: 3.69023e-09
Point: 25000 25000 3.481e+06 becomes: 0.785398 1.56076 -2.87557e+06. Back transformed: 25000 25000 3.481e+06. Error: 5.73998e-11
Point: 25000 25000 1e+06 becomes: 0.785398 1.53691 -5.35615e+06. Back transformed: 25000 25000 1e+06. Error: 9.74229e-08
Point: 25000 25000 500000 becomes: 0.785398 1.50575 -5.8556e+06. Back transformed: 25000 25000 500002. Error: 3.4858e-06
Point: 25000 25000 50000 becomes: 0.785398 1.21599 -6.29522e+06. Back transformed: 25000 25000 55272.2. Error: 0.0860949
Test for points outside of the provided depth range
Point: 6.371e+06 0 0 becomes: 0 0 -7137. Back transformed: 6.371e+06 0 0. Error: 0
Point: 6.171e+06 0 0 becomes: 0 0 -207137. Back transformed: 6.171e+06 0 0. Error: 0
Point: 3e+06 3e+06 0 becomes: 0.785398 0 -2.1355e+06. Back transformed: 3e+06 3e+06 0. Error: 1.09757e-16
Point: 3e+06 3e+06 3e+06 becomes: 0.785398 0.619368 -1.17483e+06. Back transformed: 3e+06 3e+06 3e+06. Error: 3.69023e-09
Point: 3e+06 -3e+06 -3e+06 becomes: -0.785398 -0.619368 -1.17483e+06. Back transformed: 3e+06 -3e+06 -3e+06. Error: 3.69023e-09
Point: 25000 25000 3.481e+06 becomes: 0.785398 1.56076 -2.87557e+06. Back transformed: 25000 25000 3.481e+06. Error: 5.73998e-11
Point: 25000 25000 1e+06 becomes: 0.785398 1.53691 -5.35615e+06. Back transformed: 25000 25000 1e+06. Error: 9.74229e-08
Point: 25000 25000 500000 becomes: 0.785398 1.50575 -5.8556e+06. Back transformed: 25000 25000 500002. Error: 3.4858e-06
Point: 25000 25000 50000 becomes: 0.785398 1.21599 -6.29522e+06. Back transformed: 25000 25000 55272.2. Error: 0.0860949

0 comments on commit 5a8e038

Please sign in to comment.