Line data Source code
1 : /********************************************************************/
2 : /* SOFTWARE COPYRIGHT NOTIFICATION */
3 : /* Cardinal */
4 : /* */
5 : /* (c) 2021 UChicago Argonne, LLC */
6 : /* ALL RIGHTS RESERVED */
7 : /* */
8 : /* Prepared by UChicago Argonne, LLC */
9 : /* Under Contract No. DE-AC02-06CH11357 */
10 : /* With the U. S. Department of Energy */
11 : /* */
12 : /* Prepared by Battelle Energy Alliance, LLC */
13 : /* Under Contract No. DE-AC07-05ID14517 */
14 : /* With the U. S. Department of Energy */
15 : /* */
16 : /* See LICENSE for full restrictions */
17 : /********************************************************************/
18 :
19 : #ifdef ENABLE_NEK_COUPLING
20 :
21 : #include "NekBinnedSideIntegral.h"
22 : #include "NekInterface.h"
23 :
24 : registerMooseObject("CardinalApp", NekBinnedSideIntegral);
25 :
26 : InputParameters
27 113 : NekBinnedSideIntegral::validParams()
28 : {
29 113 : InputParameters params = NekSideSpatialBinUserObject::validParams();
30 113 : params.addClassDescription(
31 : "Compute the spatially-binned side integral of a field over a boundary of the NekRS mesh");
32 113 : return params;
33 0 : }
34 :
35 57 : NekBinnedSideIntegral::NekBinnedSideIntegral(const InputParameters & parameters)
36 57 : : NekSideSpatialBinUserObject(parameters)
37 : {
38 56 : if (_fixed_mesh)
39 56 : computeBinVolumes();
40 56 : }
41 :
42 : void
43 56 : NekBinnedSideIntegral::getBinVolumes()
44 : {
45 56 : resetPartialStorage();
46 :
47 56 : mesh_t * mesh = nekrs::entireMesh();
48 56 : const auto & sgeo = nekrs::getSgeo();
49 :
50 92688 : for (int i = 0; i < mesh->Nelements; ++i)
51 : {
52 648424 : for (int j = 0; j < mesh->Nfaces; ++j)
53 : {
54 555792 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
55 555792 : if (std::find(_boundary.begin(), _boundary.end(), face_id) != _boundary.end())
56 : {
57 38472 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
58 625224 : for (int v = 0; v < mesh->Nfp; ++v)
59 : {
60 586752 : Point p = nekPoint(i, j, v);
61 586752 : unsigned int b = bin(p);
62 586752 : _bin_partial_values[b] += sgeo[mesh->Nsgeo * (offset + v) + WSJID];
63 586752 : _bin_partial_counts[b]++;
64 : }
65 : }
66 : }
67 : }
68 :
69 : // sum across all processes
70 56 : MPI_Allreduce(
71 : _bin_partial_values, _bin_volumes, _n_bins, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
72 56 : MPI_Allreduce(
73 : _bin_partial_counts, _bin_counts, _n_bins, MPI_INT, MPI_SUM, platform->comm.mpiComm);
74 :
75 : // dimensionalize
76 272 : for (unsigned int i = 0; i < _n_bins; ++i)
77 216 : nekrs::dimensionalizeArea(_bin_volumes[i]);
78 56 : }
79 :
80 : void
81 72 : NekBinnedSideIntegral::binnedSideIntegral(const field::NekFieldEnum & integrand,
82 : double * total_integral)
83 : {
84 72 : resetPartialStorage();
85 :
86 72 : mesh_t * mesh = nekrs::entireMesh();
87 72 : double (*f)(int, int) = nekrs::solutionPointer(integrand);
88 72 : const auto & sgeo = nekrs::getSgeo();
89 :
90 107936 : for (int i = 0; i < mesh->Nelements; ++i)
91 : {
92 755048 : for (int j = 0; j < mesh->Nfaces; ++j)
93 : {
94 647184 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
95 647184 : if (std::find(_boundary.begin(), _boundary.end(), face_id) != _boundary.end())
96 : {
97 39240 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
98 638280 : for (int v = 0; v < mesh->Nfp; ++v)
99 : {
100 599040 : Point p = nekPoint(i, j, v);
101 599040 : unsigned int b = bin(p);
102 599040 : _bin_partial_values[b] +=
103 599040 : f(mesh->vmapM[offset + v], 0) * sgeo[mesh->Nsgeo * (offset + v) + WSJID];
104 : }
105 : }
106 : }
107 : }
108 :
109 : // sum across all processes
110 72 : MPI_Allreduce(
111 : _bin_partial_values, total_integral, _n_bins, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
112 :
113 : // dimensionalize
114 368 : for (unsigned int i = 0; i < _n_bins; ++i)
115 296 : nekrs::dimensionalizeSideIntegral(integrand, _bin_volumes[i], total_integral[i]);
116 72 : }
117 :
118 : Real
119 0 : NekBinnedSideIntegral::spatialValue(const Point & p, const unsigned int & component) const
120 : {
121 0 : const auto & i = bin(p);
122 0 : return _bin_values[i] * _velocity_bin_directions[i](component);
123 : }
124 :
125 : void
126 72 : NekBinnedSideIntegral::computeIntegral()
127 : {
128 : // if the mesh is changing, re-compute the areas of the bins
129 72 : if (!_fixed_mesh)
130 0 : computeBinVolumes();
131 :
132 72 : if (_field == field::velocity_component)
133 : {
134 0 : binnedSideIntegral(field::velocity_x, _bin_values_x);
135 0 : binnedSideIntegral(field::velocity_y, _bin_values_y);
136 0 : binnedSideIntegral(field::velocity_z, _bin_values_z);
137 :
138 0 : for (unsigned int i = 0; i < num_bins(); ++i)
139 : {
140 0 : Point velocity(_bin_values_x[i], _bin_values_y[i], _bin_values_z[i]);
141 0 : _bin_values[i] = _velocity_bin_directions[i] * velocity;
142 : }
143 : }
144 : else
145 72 : binnedSideIntegral(_field, _bin_values);
146 72 : }
147 :
148 : void
149 72 : NekBinnedSideIntegral::executeUserObject()
150 : {
151 72 : computeIntegral();
152 72 : }
153 :
154 : #endif
|
