From 36586c0b078cefda1931665230b44395881123eb Mon Sep 17 00:00:00 2001 From: "Professor E. G. Puckett" Date: Thu, 15 Dec 2016 12:26:40 -0800 Subject: [PATCH] Added particle interpolation scheme bilinear least squares --- doc/modules/changes/20170513_egp-cig-reu | 7 + .../interpolator/bilinear_least_squares.h | 62 ++++++ .../interpolator/bilinear_least_squares.cc | 180 ++++++++++++++++++ tests/particle_interpolator_bilinear_solkz.cc | 2 + .../particle_interpolator_bilinear_solkz.prm | 116 +++++++++++ .../screen-output | 24 +++ .../statistics | 14 ++ 7 files changed, 405 insertions(+) create mode 100644 doc/modules/changes/20170513_egp-cig-reu create mode 100644 include/aspect/particle/interpolator/bilinear_least_squares.h create mode 100644 source/particle/interpolator/bilinear_least_squares.cc create mode 100644 tests/particle_interpolator_bilinear_solkz.cc create mode 100644 tests/particle_interpolator_bilinear_solkz.prm create mode 100644 tests/particle_interpolator_bilinear_solkz/screen-output create mode 100644 tests/particle_interpolator_bilinear_solkz/statistics diff --git a/doc/modules/changes/20170513_egp-cig-reu b/doc/modules/changes/20170513_egp-cig-reu new file mode 100644 index 00000000000..3504c952bdc --- /dev/null +++ b/doc/modules/changes/20170513_egp-cig-reu @@ -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. +
+(Elbridge Gerry Puckett, Ying He, Harsha Lokavarapu 2017/05/13) \ No newline at end of file diff --git a/include/aspect/particle/interpolator/bilinear_least_squares.h b/include/aspect/particle/interpolator/bilinear_least_squares.h new file mode 100644 index 00000000000..a8b15083abf --- /dev/null +++ b/include/aspect/particle/interpolator/bilinear_least_squares.h @@ -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 + . + */ + +#ifndef _aspect_particle_interpolator_bilinear_least_squares_h +#define _aspect_particle_interpolator_bilinear_least_squares_h + +#include +#include + +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 + class BilinearLeastSquares : public Interface, public aspect::SimulatorAccess + { + 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 > + properties_at_points(const std::multimap > &particles, + const std::vector > &positions, + const ComponentMask &selected_properties, + const typename parallel::distributed::Triangulation::active_cell_iterator &cell) const; + + // avoid -Woverloaded-virtual: + using Interface::properties_at_points; + }; + } + } +} + +#endif diff --git a/source/particle/interpolator/bilinear_least_squares.cc b/source/particle/interpolator/bilinear_least_squares.cc new file mode 100644 index 00000000000..b3e3c3955ed --- /dev/null +++ b/source/particle/interpolator/bilinear_least_squares.cc @@ -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 + . + */ + +#include + +#include +#include +#include + +namespace aspect +{ + namespace Particle + { + namespace Interpolator + { + template + std::vector > + BilinearLeastSquares::properties_at_points(const std::multimap > &particles, + const std::vector > &positions, + const ComponentMask &selected_properties, + const typename parallel::distributed::Triangulation::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 approximated_cell_midpoint = std::accumulate (positions.begin(), positions.end(), Point()) + / static_cast (positions.size()); + + typename parallel::distributed::Triangulation::active_cell_iterator found_cell; + + if (cell == typename parallel::distributed::Triangulation::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 (found_cell->level(),found_cell->index()); + const std::pair >::const_iterator, + typename std::multimap >::const_iterator> particle_range = particles.equal_range(cell_index); + + + std::vector > cell_properties(positions.size(), + std::vector(n_particle_properties, + numbers::signaling_nan())); + + 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 A(n_particles, matrix_dimension); + Vector 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 >::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 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 B(matrix_dimension, matrix_dimension); + dealii::LAPACKFullMatrix B_inverse(B); + + Vector c_ATr(matrix_dimension); + Vector 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 >::const_iterator itr = positions.begin(); itr != positions.end(); ++itr, ++index_positions) + { + const Point 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.") + } + } +} diff --git a/tests/particle_interpolator_bilinear_solkz.cc b/tests/particle_interpolator_bilinear_solkz.cc new file mode 100644 index 00000000000..3dd35c06451 --- /dev/null +++ b/tests/particle_interpolator_bilinear_solkz.cc @@ -0,0 +1,2 @@ +#include "../benchmarks/solkz/compositional_fields/solkz_compositional_fields.cc" + diff --git a/tests/particle_interpolator_bilinear_solkz.prm b/tests/particle_interpolator_bilinear_solkz.prm new file mode 100644 index 00000000000..3e62de9aba2 --- /dev/null +++ b/tests/particle_interpolator_bilinear_solkz.prm @@ -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 \ No newline at end of file diff --git a/tests/particle_interpolator_bilinear_solkz/screen-output b/tests/particle_interpolator_bilinear_solkz/screen-output new file mode 100644 index 00000000000..c69c4642fbf --- /dev/null +++ b/tests/particle_interpolator_bilinear_solkz/screen-output @@ -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 + + ++---------------------------------------------+------------+------------+ ++---------------------------------+-----------+------------+------------+ ++---------------------------------+-----------+------------+------------+ + diff --git a/tests/particle_interpolator_bilinear_solkz/statistics b/tests/particle_interpolator_bilinear_solkz/statistics new file mode 100644 index 00000000000..7c5077eba42 --- /dev/null +++ b/tests/particle_interpolator_bilinear_solkz/statistics @@ -0,0 +1,14 @@ +# 1: Time step number +# 2: Time (years) +# 3: Time step size (years) +# 4: Number of mesh cells +# 5: Number of Stokes degrees of freedom +# 6: Number of temperature degrees of freedom +# 7: Number of degrees of freedom for all compositions +# 8: Iterations for temperature solver +# 9: Iterations for Stokes solver +# 10: Velocity iterations in Stokes preconditioner +# 11: Schur complement iterations in Stokes preconditioner +# 12: Number of advected particles +# 13: Particle file name +0 0.000000000000e+00 0.000000000000e+00 1024 9539 4225 18432 0 46 47 46 16384 output-particle_interpolator_bilinear_solkz/output-particle_interpolator_bilinear_solkz/particles/particles-00000