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 "NekBinnedVolumeIntegral.h" 22 : #include "CardinalUtils.h" 23 : #include "NekInterface.h" 24 : 25 : registerMooseObject("CardinalApp", NekBinnedVolumeIntegral); 26 : 27 : InputParameters 28 261 : NekBinnedVolumeIntegral::validParams() 29 : { 30 261 : InputParameters params = NekVolumeSpatialBinUserObject::validParams(); 31 261 : params.addClassDescription( 32 : "Compute the spatially-binned volume integral of a field over the NekRS mesh"); 33 261 : return params; 34 0 : } 35 : 36 133 : NekBinnedVolumeIntegral::NekBinnedVolumeIntegral(const InputParameters & parameters) 37 133 : : NekVolumeSpatialBinUserObject(parameters) 38 : { 39 126 : if (_fixed_mesh) 40 126 : computeBinVolumes(); 41 125 : } 42 : 43 : void 44 126 : NekBinnedVolumeIntegral::getBinVolumes() 45 : { 46 126 : resetPartialStorage(); 47 : 48 126 : mesh_t * mesh = nekrs::entireMesh(); 49 126 : const auto & vgeo = nekrs::getVgeo(); 50 141414 : for (int k = 0; k < mesh->Nelements; ++k) 51 : { 52 141288 : int offset = k * mesh->Np; 53 13137768 : for (int v = 0; v < mesh->Np; ++v) 54 : { 55 12996480 : Point p = nekPoint(k, v); 56 12996480 : unsigned int b = bin(p); 57 12996480 : _bin_partial_values[b] += vgeo[mesh->Nvgeo * offset + v + mesh->Np * JWID]; 58 12996480 : _bin_partial_counts[b]++; 59 : } 60 : } 61 : 62 : // sum across all processes 63 126 : MPI_Allreduce( 64 : _bin_partial_values, _bin_volumes, _n_bins, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm); 65 126 : MPI_Allreduce( 66 : _bin_partial_counts, _bin_counts, _n_bins, MPI_INT, MPI_SUM, platform->comm.mpiComm); 67 : 68 : // dimensionalize 69 17049 : for (unsigned int i = 0; i < _n_bins; ++i) 70 16923 : nekrs::dimensionalizeVolume(_bin_volumes[i]); 71 126 : } 72 : 73 : Real 74 20952 : NekBinnedVolumeIntegral::spatialValue(const Point & p, const unsigned int & component) const 75 : { 76 20952 : const auto & i = bin(p); 77 20952 : return _bin_values[i] * _velocity_bin_directions[i](component); 78 : } 79 : 80 : void 81 140 : NekBinnedVolumeIntegral::binnedVolumeIntegral(const field::NekFieldEnum & integrand, 82 : double * total_integral) 83 : { 84 140 : resetPartialStorage(); 85 : 86 140 : mesh_t * mesh = nekrs::entireMesh(); 87 140 : double (*f)(int, int) = nekrs::solutionPointer(integrand); 88 140 : const auto & vgeo = nekrs::getVgeo(); 89 : 90 155612 : for (int k = 0; k < mesh->Nelements; ++k) 91 : { 92 155472 : int offset = k * mesh->Np; 93 14156496 : for (int v = 0; v < mesh->Np; ++v) 94 : { 95 14001024 : Point p = nekPoint(k, v); 96 14001024 : unsigned int b = bin(p); 97 14001024 : _bin_partial_values[b] += f(offset + v, 0) * vgeo[mesh->Nvgeo * offset + v + mesh->Np * JWID]; 98 : } 99 : } 100 : 101 : // sum across all processes 102 140 : MPI_Allreduce( 103 : _bin_partial_values, total_integral, _n_bins, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm); 104 : 105 8948 : for (unsigned int i = 0; i < _n_bins; ++i) 106 8808 : nekrs::dimensionalizeVolumeIntegral(integrand, _bin_volumes[i], total_integral[i]); 107 140 : } 108 : 109 : void 110 132 : NekBinnedVolumeIntegral::executeUserObject() 111 : { 112 : // if the mesh is changing, re-compute the volumes of the bins and check the counts 113 132 : if (!_fixed_mesh) 114 0 : computeBinVolumes(); 115 : 116 132 : if (_field == field::velocity_component) 117 : { 118 4 : binnedVolumeIntegral(field::velocity_x, _bin_values_x); 119 4 : binnedVolumeIntegral(field::velocity_y, _bin_values_y); 120 4 : binnedVolumeIntegral(field::velocity_z, _bin_values_z); 121 : 122 436 : for (unsigned int i = 0; i < num_bins(); ++i) 123 : { 124 432 : Point velocity(_bin_values_x[i], _bin_values_y[i], _bin_values_z[i]); 125 432 : _bin_values[i] = _velocity_bin_directions[i] * velocity; 126 : } 127 : } 128 : else 129 128 : binnedVolumeIntegral(_field, _bin_values); 130 132 : } 131 : 132 : #endif