Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #1554 from hlokavarapu/particle_interpolation_sche…
…me_bilinear Bilinear Least Squares Particle Interpolation Scheme
- Loading branch information
Showing
7 changed files
with
405 additions
and
0 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,7 @@ | ||
New particle interpolation scheme bilinear least squares | ||
using singular value decomposition algorithm. Currently | ||
only 2D models are supported. We chose a simple overshoot | ||
and undershoot correction scheme, mainly to truncate based | ||
on local cell maximum and minimum property value. | ||
<br> | ||
(Elbridge Gerry Puckett, Ying He, Harsha Lokavarapu 2017/05/13) |
62 changes: 62 additions & 0 deletions
62
include/aspect/particle/interpolator/bilinear_least_squares.h
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,62 @@ | ||
/* | ||
Copyright (C) 2017 by the authors of the ASPECT code. | ||
This file is part of ASPECT. | ||
ASPECT is free software; you can redistribute it and/or modify | ||
it under the terms of the GNU General Public License as published by | ||
the Free Software Foundation; either version 2, or (at your option) | ||
any later version. | ||
ASPECT is distributed in the hope that it will be useful, | ||
but WITHOUT ANY WARRANTY; without even the implied warranty of | ||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | ||
GNU General Public License for more details. | ||
You should have received a copy of the GNU General Public License | ||
along with ASPECT; see the file doc/COPYING. If not see | ||
<http://www.gnu.org/licenses/>. | ||
*/ | ||
|
||
#ifndef _aspect_particle_interpolator_bilinear_least_squares_h | ||
#define _aspect_particle_interpolator_bilinear_least_squares_h | ||
|
||
#include <aspect/particle/interpolator/interface.h> | ||
#include <aspect/simulator_access.h> | ||
|
||
namespace aspect | ||
{ | ||
namespace Particle | ||
{ | ||
namespace Interpolator | ||
{ | ||
/** | ||
* Return the interpolated properties of all particles of the given cell using bilinear least squares method. | ||
* Currently, only the two dimensional model is supported. | ||
* | ||
* @ingroup ParticleInterpolators | ||
*/ | ||
template <int dim> | ||
class BilinearLeastSquares : public Interface<dim>, public aspect::SimulatorAccess<dim> | ||
{ | ||
public: | ||
/** | ||
* Return the cell-wise evaluated properties of the bilinear least squares function at the positions. | ||
* | ||
* @copydoc aspect::Particle::Interpolator::Interface::properties_at_points() | ||
*/ | ||
virtual | ||
std::vector<std::vector<double> > | ||
properties_at_points(const std::multimap<types::LevelInd, Particle<dim> > &particles, | ||
const std::vector<Point<dim> > &positions, | ||
const ComponentMask &selected_properties, | ||
const typename parallel::distributed::Triangulation<dim>::active_cell_iterator &cell) const; | ||
|
||
// avoid -Woverloaded-virtual: | ||
using Interface<dim>::properties_at_points; | ||
}; | ||
} | ||
} | ||
} | ||
|
||
#endif |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,180 @@ | ||
/* | ||
Copyright (C) 2017 by the authors of the ASPECT code. | ||
This file is part of ASPECT. | ||
ASPECT is free software; you can redistribute it and/or modify | ||
it under the terms of the GNU General Public License as published by | ||
the Free Software Foundation; either version 2, or (at your option) | ||
any later version. | ||
ASPECT is distributed in the hope that it will be useful, | ||
but WITHOUT ANY WARRANTY; without even the implied warranty of | ||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | ||
GNU General Public License for more details. | ||
You should have received a copy of the GNU General Public License | ||
along with ASPECT; see the file doc/COPYING. If not see | ||
<http://www.gnu.org/licenses/>. | ||
*/ | ||
|
||
#include <aspect/particle/interpolator/bilinear_least_squares.h> | ||
|
||
#include <deal.II/grid/grid_tools.h> | ||
#include <deal.II/base/signaling_nan.h> | ||
#include <deal.II/lac/full_matrix.templates.h> | ||
|
||
namespace aspect | ||
{ | ||
namespace Particle | ||
{ | ||
namespace Interpolator | ||
{ | ||
template <int dim> | ||
std::vector<std::vector<double> > | ||
BilinearLeastSquares<dim>::properties_at_points(const std::multimap<types::LevelInd, Particle<dim> > &particles, | ||
const std::vector<Point<dim> > &positions, | ||
const ComponentMask &selected_properties, | ||
const typename parallel::distributed::Triangulation<dim>::active_cell_iterator &cell) const | ||
{ | ||
const unsigned int n_particle_properties = particles.begin()->second.get_properties().size(); | ||
|
||
const unsigned int property_index = selected_properties.first_selected_component(selected_properties.size()); | ||
|
||
AssertThrow(property_index != numbers::invalid_unsigned_int, | ||
ExcMessage("Internal error: the particle property interpolator was " | ||
"called without a specified component to interpolate.")); | ||
|
||
AssertThrow(dim == 2, | ||
ExcMessage("Currently, the particle interpolator 'bilinear' is only supported for 2D models.")); | ||
|
||
AssertThrow(selected_properties.n_selected_components(n_particle_properties) == 1, | ||
ExcNotImplemented("Interpolation of multiple components is not supported.")); | ||
|
||
const Point<dim> approximated_cell_midpoint = std::accumulate (positions.begin(), positions.end(), Point<dim>()) | ||
/ static_cast<double> (positions.size()); | ||
|
||
typename parallel::distributed::Triangulation<dim>::active_cell_iterator found_cell; | ||
|
||
if (cell == typename parallel::distributed::Triangulation<dim>::active_cell_iterator()) | ||
{ | ||
// We can not simply use one of the points as input for find_active_cell_around_point | ||
// because for vertices of mesh cells we might end up getting ghost_cells as return value | ||
// instead of the local active cell. So make sure we are well in the inside of a cell. | ||
Assert(positions.size() > 0, | ||
ExcMessage("The particle property interpolator was not given any " | ||
"positions to evaluate the particle cell_properties at.")); | ||
|
||
|
||
found_cell = | ||
(GridTools::find_active_cell_around_point<> (this->get_mapping(), | ||
this->get_triangulation(), | ||
approximated_cell_midpoint)).first; | ||
} | ||
else | ||
found_cell = cell; | ||
|
||
const types::LevelInd cell_index = std::make_pair<unsigned int, unsigned int> (found_cell->level(),found_cell->index()); | ||
const std::pair<typename std::multimap<types::LevelInd, Particle<dim> >::const_iterator, | ||
typename std::multimap<types::LevelInd, Particle<dim> >::const_iterator> particle_range = particles.equal_range(cell_index); | ||
|
||
|
||
std::vector<std::vector<double> > cell_properties(positions.size(), | ||
std::vector<double>(n_particle_properties, | ||
numbers::signaling_nan<double>())); | ||
|
||
const unsigned int n_particles = std::distance(particle_range.first,particle_range.second); | ||
|
||
AssertThrow(n_particles != 0, | ||
ExcMessage("At least one cell contained no particles. The 'bilinear'" | ||
"interpolation scheme does not support this case. ")); | ||
|
||
|
||
// Noticed that the size of matrix A is n_particles x matrix_dimension | ||
// which usually is not a square matrix. Therefore, we solve Ax=r by | ||
// solving A^TAx= A^Tr. | ||
const unsigned int matrix_dimension = 4; | ||
dealii::LAPACKFullMatrix<double> A(n_particles, matrix_dimension); | ||
Vector<double> r(n_particles); | ||
r = 0; | ||
|
||
// Keep track of the local min and max in order to correct | ||
// overshoot and undershoot of the interpolated particle property. | ||
double max_value_for_particle_property = (particle_range.first)->second.get_properties()[property_index]; | ||
double min_value_for_particle_property = (particle_range.first)->second.get_properties()[property_index]; | ||
|
||
unsigned int index = 0; | ||
const double cell_diameter = found_cell->diameter(); | ||
for (typename std::multimap<types::LevelInd, Particle<dim> >::const_iterator particle = particle_range.first; | ||
particle != particle_range.second; ++particle, ++index) | ||
{ | ||
const double particle_property_value = particle->second.get_properties()[property_index]; | ||
r[index] = particle_property_value; | ||
|
||
const Point<dim> position = particle->second.get_location(); | ||
A(index,0) = 1; | ||
A(index,1) = (position[0] - approximated_cell_midpoint[0])/cell_diameter; | ||
A(index,2) = (position[1] - approximated_cell_midpoint[1])/cell_diameter; | ||
A(index,3) = (position[0] - approximated_cell_midpoint[0]) * (position[1] - approximated_cell_midpoint[1])/std::pow(cell_diameter,2); | ||
|
||
if (min_value_for_particle_property > particle_property_value) | ||
min_value_for_particle_property = particle_property_value; | ||
if (max_value_for_particle_property < particle_property_value) | ||
max_value_for_particle_property = particle_property_value; | ||
} | ||
|
||
dealii::LAPACKFullMatrix<double> B(matrix_dimension, matrix_dimension); | ||
dealii::LAPACKFullMatrix<double> B_inverse(B); | ||
|
||
Vector<double> c_ATr(matrix_dimension); | ||
Vector<double> c(matrix_dimension); | ||
|
||
const double threshold = 1e-15; | ||
unsigned int index_positions = 0; | ||
|
||
// Matrix A can be rank deficient if it does not have full rank, therefore singular. | ||
// To circumvent this issue, we solve A^TAx=A^Tr by using singular value | ||
// decomposition (SVD). | ||
A.Tmmult(B, A, false); | ||
A.Tvmult(c_ATr,r); | ||
B_inverse.compute_inverse_svd(threshold); | ||
B_inverse.vmult(c, c_ATr); | ||
|
||
for (typename std::vector<Point<dim> >::const_iterator itr = positions.begin(); itr != positions.end(); ++itr, ++index_positions) | ||
{ | ||
const Point<dim> support_point = *itr; | ||
double interpolated_value = c[0] + | ||
c[1]*(support_point[0] - approximated_cell_midpoint[0])/cell_diameter + | ||
c[2]*(support_point[1] - approximated_cell_midpoint[1])/cell_diameter + | ||
c[3]*(support_point[0] - approximated_cell_midpoint[0])*(support_point[1] - approximated_cell_midpoint[1])/std::pow(cell_diameter,2); | ||
|
||
cell_properties[index_positions][property_index] = interpolated_value; | ||
|
||
// Overshoot and undershoot correction of interpolated particle property. | ||
if (interpolated_value > max_value_for_particle_property) | ||
interpolated_value = max_value_for_particle_property; | ||
else if (interpolated_value < min_value_for_particle_property) | ||
interpolated_value = min_value_for_particle_property; | ||
} | ||
return cell_properties; | ||
} | ||
} | ||
} | ||
} | ||
|
||
|
||
// explicit instantiations | ||
namespace aspect | ||
{ | ||
namespace Particle | ||
{ | ||
namespace Interpolator | ||
{ | ||
ASPECT_REGISTER_PARTICLE_INTERPOLATOR(BilinearLeastSquares, | ||
"bilinear least squares", | ||
"Interpolates particle properties onto a vector of points using a " | ||
"bilinear least squares method. Currently only 2D models are " | ||
"supported. Note that deal.II must be configured with BLAS/LAPACK.") | ||
} | ||
} | ||
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,2 @@ | ||
#include "../benchmarks/solkz/compositional_fields/solkz_compositional_fields.cc" | ||
|
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,116 @@ | ||
# A test based on the SolKz benchmark that performs bilinear | ||
# least squares interpolation from particle properties to | ||
# compositional fields that are flagged as 'particle' advected fields. | ||
|
||
set Additional shared libraries = ./libparticle_interpolator_bilinear_solkz.so | ||
|
||
############### Global parameters | ||
|
||
set Dimension = 2 | ||
|
||
set Start time = 0 | ||
set End time = 0 | ||
|
||
set Output directory = output | ||
|
||
set Pressure normalization = volume | ||
|
||
############### Parameters describing the model | ||
|
||
subsection Geometry model | ||
set Model name = box | ||
|
||
subsection Box | ||
set X extent = 1 | ||
set Y extent = 1 | ||
end | ||
end | ||
|
||
subsection Model settings | ||
set Prescribed velocity boundary indicators = | ||
set Tangential velocity boundary indicators = left, right, bottom, top | ||
set Zero velocity boundary indicators = | ||
end | ||
|
||
subsection Material model | ||
set Model name = SolKzCompositionalMaterial | ||
end | ||
|
||
subsection Gravity model | ||
set Model name = vertical | ||
end | ||
|
||
############### Parameters describing the temperature field | ||
|
||
subsection Boundary temperature model | ||
set Model name = box | ||
end | ||
|
||
subsection Initial temperature model | ||
set Model name = perturbed box | ||
end | ||
|
||
############### Parameters describing the discretization | ||
|
||
subsection Discretization | ||
set Stokes velocity polynomial degree = 2 | ||
set Use locally conservative discretization = false | ||
set Use discontinuous composition discretization = true | ||
end | ||
|
||
subsection Mesh refinement | ||
set Initial adaptive refinement = 0 | ||
set Initial global refinement = 5 | ||
end | ||
|
||
############### Parameters describing the compositional field | ||
# Note: The compositional field is what drives the flow | ||
# in this example | ||
|
||
subsection Compositional fields | ||
set Number of fields = 2 | ||
set Names of fields = density_comp, viscosity_comp | ||
set Compositional field methods = particles, particles | ||
set Mapped particle properties = density_comp:function[0], viscosity_comp:function[1] | ||
end | ||
|
||
subsection Initial composition model | ||
set Model name = function | ||
subsection Function | ||
set Variable names = x,z | ||
set Function constants = pi=3.1415926, eta_b=1e6 | ||
set Function expression = -1 * sin(2*z) * cos(3*pi*x); exp(log(eta_b)*z) | ||
end | ||
end | ||
|
||
############### Parameters describing what to do with the solution | ||
|
||
subsection Postprocess | ||
set List of postprocessors = particles, SolKzPostprocessor | ||
|
||
subsection Particles | ||
set Number of particles = 16384 | ||
set Time between data output = 0 | ||
set Data output format = vtu | ||
set List of particle properties = function | ||
set Integration scheme = rk2 | ||
set Interpolation scheme = cell average | ||
set Maximum particles per cell = 16384 | ||
set Update ghost particles = true | ||
|
||
subsection Function | ||
set Number of components = 2 | ||
set Variable names = x, z | ||
set Function constants = pi=3.1415926, eta_b=1e6 | ||
set Function expression = -1 * sin(2*z) * cos(3*pi*x); exp(log(eta_b)*z) | ||
end | ||
|
||
set Particle generator name = reference cell | ||
|
||
subsection Generator | ||
subsection Reference cell | ||
set Number of particles per cell per direction = 4 | ||
end | ||
end | ||
end | ||
end |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,24 @@ | ||
----------------------------------------------------------------------------- | ||
----------------------------------------------------------------------------- | ||
|
||
Loading shared library <./libparticle_interpolator_bilinear_solkz.so> | ||
|
||
Number of active cells: 1,024 (on 6 levels) | ||
Number of degrees of freedom: 32,196 (8,450+1,089+4,225+9,216+9,216) | ||
|
||
*** Timestep 0: t=0 years | ||
Skipping temperature solve because RHS is zero. | ||
Rebuilding Stokes preconditioner... | ||
Solving Stokes system... 46+0 iterations. | ||
|
||
Postprocessing: | ||
Writing particle output: output-particle_interpolator_bilinear_solkz/particles/particles-00000 | ||
Errors u_L1, p_L1, u_L2, p_L2: 2.396693e-07, 3.813761e-04, 3.798509e-07, 1.430329e-03 | ||
|
||
Termination requested by criterion: end time | ||
|
||
|
||
+---------------------------------------------+------------+------------+ | ||
+---------------------------------+-----------+------------+------------+ | ||
+---------------------------------+-----------+------------+------------+ | ||
|
Oops, something went wrong.