Skip to content

Commit 792e91f

Browse files
gassmoellergassmoeller
committed
Add distance weighted average interpolator
1 parent 0d485a6 commit 792e91f

6 files changed

Lines changed: 398 additions & 0 deletions

File tree

‎doc/modules/changes/20240609_gassmoeller‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,8 @@
1+
Added: A new particle interpolator plugin 'distance weighted average' that
2+
computes the interpolated properties as an average of particle properties
3+
weighted by a distance function. This plugin is similarly accurate
4+
as the 'cell average' interpolator, but provides smoother results and is
5+
less susceptible to jumps in interpolated properties if particles move
6+
small distances.
7+
<br>
8+
(Rene Gassmoeller, 2024/06/09)

‎include/aspect/particle/interpolator/distance_weighted_average.h‎

Lines changed: 72 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,72 @@
1+
/*
2+
Copyright (C) 2024 by the authors of the ASPECT code.
3+
4+
This file is part of ASPECT.
5+
6+
ASPECT is free software; you can redistribute it and/or modify
7+
it under the terms of the GNU General Public License as published by
8+
the Free Software Foundation; either version 2, or (at your option)
9+
any later version.
10+
11+
ASPECT is distributed in the hope that it will be useful,
12+
but WITHOUT ANY WARRANTY; without even the implied warranty of
13+
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14+
GNU General Public License for more details.
15+
16+
You should have received a copy of the GNU General Public License
17+
along with ASPECT; see the file LICENSE. If not see
18+
<http://www.gnu.org/licenses/>.
19+
*/
20+
21+
#ifndef _aspect_particle_interpolator_distance_weighted_average_h
22+
#define _aspect_particle_interpolator_distance_weighted_average_h
23+
24+
#include <aspect/particle/interpolator/interface.h>
25+
#include <aspect/simulator_access.h>
26+
27+
#include <deal.II/grid/grid_tools_cache.h>
28+
29+
namespace aspect
30+
{
31+
namespace Particle
32+
{
33+
namespace Interpolator
34+
{
35+
/**
36+
* Return the interpolated properties of all particles of the given cell using bilinear least squares method.
37+
* Currently, only the two dimensional model is supported.
38+
*
39+
* @ingroup ParticleInterpolators
40+
*/
41+
template <int dim>
42+
class DistanceWeightedAverage : public Interface<dim>, public aspect::SimulatorAccess<dim>
43+
{
44+
public:
45+
/**
46+
* Update function. Called once at the beginning of each time step.
47+
*/
48+
void update() override;
49+
50+
/**
51+
* Return the cell-wise evaluated particle properties at the given @p positions.
52+
*/
53+
std::vector<std::vector<double>>
54+
properties_at_points(const ParticleHandler<dim> &particle_handler,
55+
const std::vector<Point<dim>> &positions,
56+
const ComponentMask &selected_properties,
57+
const typename parallel::distributed::Triangulation<dim>::active_cell_iterator &cell) const override;
58+
59+
// avoid -Woverloaded-virtual:
60+
using Interface<dim>::properties_at_points;
61+
62+
private:
63+
/**
64+
* Cached information.
65+
*/
66+
std::vector<std::set<typename Triangulation<dim>::active_cell_iterator>> vertex_to_cell_map;
67+
};
68+
}
69+
}
70+
}
71+
72+
#endif

‎source/particle/interpolator/distance_weighted_average.cc‎

Lines changed: 139 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,139 @@
1+
/*
2+
Copyright (C) 2024 by the authors of the ASPECT code.
3+
4+
This file is part of ASPECT.
5+
6+
ASPECT is free software; you can redistribute it and/or modify
7+
it under the terms of the GNU General Public License as published by
8+
the Free Software Foundation; either version 2, or (at your option)
9+
any later version.
10+
11+
ASPECT is distributed in the hope that it will be useful,
12+
but WITHOUT ANY WARRANTY; without even the implied warranty of
13+
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14+
GNU General Public License for more details.
15+
16+
You should have received a copy of the GNU General Public License
17+
along with ASPECT; see the file LICENSE. If not see
18+
<http://www.gnu.org/licenses/>.
19+
*/
20+
21+
#include <aspect/particle/interpolator/distance_weighted_average.h>
22+
23+
#include <deal.II/base/signaling_nan.h>
24+
#include <deal.II/grid/grid_tools.h>
25+
26+
namespace aspect
27+
{
28+
namespace Particle
29+
{
30+
namespace Interpolator
31+
{
32+
namespace internal
33+
{
34+
double weight (const double distance, const double cell_width)
35+
{
36+
// linear weight function (hat function)
37+
return std::max(1.0 - (distance/cell_width),0.0);
38+
}
39+
}
40+
41+
42+
43+
template <int dim>
44+
void
45+
DistanceWeightedAverage<dim>::update()
46+
{
47+
vertex_to_cell_map = GridTools::vertex_to_cell_map(this->get_triangulation());
48+
}
49+
50+
51+
52+
template <int dim>
53+
std::vector<std::vector<double>>
54+
DistanceWeightedAverage<dim>::properties_at_points(const ParticleHandler<dim> &particle_handler,
55+
const std::vector<Point<dim>> &positions,
56+
const ComponentMask &selected_properties,
57+
const typename parallel::distributed::Triangulation<dim>::active_cell_iterator &cell) const
58+
{
59+
const unsigned int n_interpolate_positions = positions.size();
60+
const unsigned int n_particle_properties = particle_handler.n_properties_per_particle();
61+
62+
// Create with signaling NaNs
63+
std::vector<std::vector<double>> cell_properties(n_interpolate_positions,
64+
std::vector<double>(n_particle_properties,
65+
numbers::signaling_nan<double>()));
66+
67+
// Set requested properties to 0.0
68+
for (unsigned int index_positions = 0; index_positions < n_interpolate_positions; ++index_positions)
69+
for (unsigned int index_properties = 0; index_properties < n_particle_properties; ++index_properties)
70+
if (selected_properties[index_properties])
71+
cell_properties[index_positions][index_properties] = 0.0;
72+
73+
std::set<typename Triangulation<dim>::active_cell_iterator> cell_and_neighbors;
74+
75+
for (const auto v : cell->vertex_indices())
76+
{
77+
const unsigned int vertex_index = cell->vertex_index(v);
78+
cell_and_neighbors.insert(vertex_to_cell_map[vertex_index].begin(),
79+
vertex_to_cell_map[vertex_index].end());
80+
}
81+
82+
const double interpolation_support = 0.5 * cell->diameter();
83+
std::vector<double> integrated_weight(n_interpolate_positions,0.0);
84+
85+
for (const auto &current_cell: cell_and_neighbors)
86+
{
87+
const typename ParticleHandler<dim>::particle_iterator_range particle_range =
88+
particle_handler.particles_in_cell(current_cell);
89+
90+
for (const auto &particle: particle_range)
91+
{
92+
const ArrayView<const double> particle_properties = particle.get_properties();
93+
unsigned int index_positions = 0;
94+
95+
for (const auto &interpolation_point: positions)
96+
{
97+
const double distance = particle.get_location().distance(interpolation_point);
98+
const double weight = internal::weight(distance,interpolation_support);
99+
100+
for (unsigned int index_properties = 0; index_properties < particle_properties.size(); ++index_properties)
101+
if (selected_properties[index_properties])
102+
cell_properties[index_positions][index_properties] += weight * particle_properties[index_properties];
103+
104+
integrated_weight[index_positions] += weight;
105+
++index_positions;
106+
}
107+
}
108+
}
109+
110+
for (unsigned int index_positions = 0; index_positions < n_interpolate_positions; ++index_positions)
111+
{
112+
AssertThrow(integrated_weight[index_positions] > 0.0, ExcInternalError());
113+
114+
for (unsigned int index_properties = 0; index_properties < n_particle_properties; ++index_properties)
115+
if (selected_properties[index_properties])
116+
cell_properties[index_positions][index_properties] /= integrated_weight[index_positions];
117+
}
118+
119+
return cell_properties;
120+
}
121+
}
122+
}
123+
}
124+
125+
126+
// explicit instantiations
127+
namespace aspect
128+
{
129+
namespace Particle
130+
{
131+
namespace Interpolator
132+
{
133+
ASPECT_REGISTER_PARTICLE_INTERPOLATOR(DistanceWeightedAverage,
134+
"distance weighted average",
135+
"Interpolates particle properties onto a vector of points using a "
136+
"distance weighed averaging method.")
137+
}
138+
}
139+
}

‎tests/particle_interpolator_distance_weighted_average.prm‎

Lines changed: 122 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,122 @@
1+
# A test that performs interpolation from particle properties to
2+
# compositional fields that are flagged as 'particle' advected
3+
# fields.
4+
5+
# MPI: 2
6+
7+
set Dimension = 2
8+
set End time = 70
9+
set Use years in output instead of seconds = false
10+
11+
subsection Geometry model
12+
set Model name = box
13+
14+
subsection Box
15+
set X extent = 0.9142
16+
set Y extent = 1.0000
17+
end
18+
end
19+
20+
subsection Boundary velocity model
21+
set Tangential velocity boundary indicators = left, right
22+
set Zero velocity boundary indicators = bottom, top
23+
end
24+
25+
subsection Material model
26+
set Model name = simple
27+
28+
subsection Simple model
29+
set Reference density = 1010
30+
set Viscosity = 1e2
31+
set Thermal expansion coefficient = 0
32+
set Density differential for compositional field 1 = -10
33+
end
34+
end
35+
36+
subsection Gravity model
37+
set Model name = vertical
38+
39+
subsection Vertical
40+
set Magnitude = 10
41+
end
42+
end
43+
44+
############### Parameters describing the temperature field
45+
# Note: The temperature plays no role in this model
46+
47+
subsection Initial temperature model
48+
set Model name = function
49+
50+
subsection Function
51+
set Function expression = 0
52+
end
53+
end
54+
55+
############### Parameters describing the compositional field
56+
# Note: The compositional field is what drives the flow
57+
# in this example
58+
59+
subsection Compositional fields
60+
set Number of fields = 3
61+
set Names of fields = advection_field, advection_particle, advection_particle2
62+
set Compositional field methods = field, particles, particles
63+
set Mapped particle properties = advection_particle2:velocity [1], advection_particle:function
64+
end
65+
66+
subsection Initial composition model
67+
set Model name = function
68+
69+
subsection Function
70+
set Variable names = x,z
71+
set Function constants = pi=3.1415926
72+
set Function expression = 0.5*(1+tanh((0.2+0.02*cos(pi*x/0.9142)-z)/0.02));0.0;0.0
73+
end
74+
end
75+
76+
subsection Discretization
77+
set Use discontinuous composition discretization = true
78+
end
79+
80+
############### Parameters describing the discretization
81+
82+
subsection Mesh refinement
83+
set Initial adaptive refinement = 0
84+
set Strategy = composition
85+
set Initial global refinement = 4
86+
set Time steps between mesh refinement = 0
87+
set Coarsening fraction = 0.05
88+
set Refinement fraction = 0.3
89+
end
90+
91+
############### Parameters describing what to do with the solution
92+
93+
subsection Postprocess
94+
set List of postprocessors = velocity statistics, composition statistics,particles
95+
96+
subsection Visualization
97+
set Interpolate output = false
98+
set Time between graphical output = 70
99+
end
100+
101+
subsection Particles
102+
set Number of particles = 100000
103+
set Time between data output = 70
104+
set Data output format = none
105+
set List of particle properties = velocity, function, initial composition
106+
set Interpolation scheme = distance weighted average
107+
set Update ghost particles = true
108+
set Particle generator name = random uniform
109+
110+
subsection Function
111+
set Variable names = x,z
112+
set Function expression = 0.5*(1+tanh((0.2+0.02*cos(pi*x/0.9142)-z)/0.02))
113+
end
114+
115+
subsection Generator
116+
subsection Probability density function
117+
set Variable names = x,z
118+
set Function expression = x*x*z
119+
end
120+
end
121+
end
122+
end

‎tests/particle_interpolator_distance_weighted_average/screen-output‎

Lines changed: 31 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,31 @@
1+
2+
Number of active cells: 256 (on 5 levels)
3+
Number of degrees of freedom: 10,468 (2,178+289+1,089+2,304+2,304+2,304)
4+
5+
*** Timestep 0: t=0 seconds, dt=0 seconds
6+
Skipping temperature solve because RHS is zero.
7+
Advecting particles... done.
8+
Solving advection_field system ... 0 iterations.
9+
Rebuilding Stokes preconditioner...
10+
Solving Stokes system... 34+0 iterations.
11+
12+
Postprocessing:
13+
RMS, max velocity: 0.000181 m/s, 0.000404 m/s
14+
Compositions min/max/mass: 0/1/0.1825 // 0/1/0.1828 // 0/0/0
15+
Number of advected particles: 100000
16+
17+
*** Timestep 1: t=70 seconds, dt=70 seconds
18+
Skipping temperature solve because RHS is zero.
19+
Advecting particles... done.
20+
Solving advection_field system ... 5 iterations.
21+
Solving Stokes system... 34+0 iterations.
22+
23+
Postprocessing:
24+
RMS, max velocity: 0.000328 m/s, 0.000732 m/s
25+
Compositions min/max/mass: -0.007595/1.023/0.1825 // 0/1/0.1829 // -0.0002466/0.0003166/-1.013e-08
26+
Number of advected particles: 100000
27+
28+
Termination requested by criterion: end time
29+
30+
31+

0 commit comments

Comments
 (0)