Line data    Source code

       1             : // This file is part of the Mamico project. For conditions of distribution
       2             : // and use, please see the copyright notice in Mamico's main folder
       3             : 
       4           0 : template <unsigned int dim, coupling::indexing::IndexTrait... scope> void coupling::filtering::NLM<dim, scope...>::operator()() {
       5             : #ifdef NLM_DEBUG
       6             :   std::cout << "    NLM: called operator()" << std::endl;
       7             : #endif
       8             : 
       9           0 :   save_data();
      10           0 :   update_patchfield();
      11           0 :   denoise();
      12           0 :   increment_time();
      13           0 : }
      14             : 
      15           0 : template <unsigned int dim, coupling::indexing::IndexTrait... scope> void coupling::filtering::NLM<dim, scope...>::save_data() {
      16           0 :   for (unsigned int i = 0; i < this->_inputCells.size(); ++i) {
      17             :     using namespace coupling::indexing;
      18           0 :     CellIndex_T idx{CellIndex<dim, IndexTrait::local, IndexTrait::md2macro, IndexTrait::noGhost>{i}};
      19           0 :     tarch::la::Vector<dim, unsigned int> idx_raw{idx.get()};
      20           0 :     coupling::filtering::Quantities<dim>& q = _flowfield(idx_raw, _t);
      21           0 :     coupling::filtering::Quantities<dim>& qp = _flowfield_prefiltered(idx_raw, _t);
      22           0 :     q[0] = this->_inputCellVectors[0][i]->getMacroscopicMass();
      23           0 :     qp[0] = this->_inputCellVectors[1][i]->getMacroscopicMass();
      24           0 :     for (unsigned int d = 0; d < dim; d++) {
      25           0 :       q[d + 1] = this->_inputCellVectors[0][i]->getMacroscopicMomentum()[d];
      26           0 :       qp[d + 1] = this->_inputCellVectors[1][i]->getMacroscopicMomentum()[d];
      27             :     }
      28             :   }
      29           0 : }
      30             : 
      31           0 : template <unsigned int dim, coupling::indexing::IndexTrait... scope> void coupling::filtering::NLM<dim, scope...>::update_patchfield() {
      32           0 :   for (auto idx : _innerCellIndices) {
      33           0 :     auto idx_inner = (tarch::la::Vector<dim, unsigned int>)(idx.get()) - tarch::la::Vector<dim, unsigned int>{1};
      34           0 :     if (_cycleCounter > _t)
      35           0 :       _patchfield.destroy(idx_inner, _t);
      36           0 :     _patchfield.construct(idx_inner, _t, tarch::la::Vector<dim, unsigned int>(2 * _d + 1), 2 * _d + 1, _flowfield_prefiltered,
      37           0 :                           (tarch::la::Vector<dim, unsigned int>)(idx.get()), _t);
      38             :   }
      39           0 : }
      40             : 
      41           0 : template <unsigned int dim, coupling::indexing::IndexTrait... scope> void coupling::filtering::NLM<dim, scope...>::denoise() {
      42             :   // counter for adaptive parameter modification
      43           0 :   int cnt_sigsq = 0, cnt_hsq = 0, cnt_total = 0;
      44           0 :   for (auto const& idx : _innerCellIndices) {
      45           0 :     auto idx_inner = (tarch::la::Vector<dim, unsigned int>)(idx.get()) - tarch::la::Vector<dim, unsigned int>{1};
      46           0 :     auto& me = _patchfield(idx_inner, _t);
      47           0 :     double C = 0;
      48             : 
      49           0 :     auto res_offset = compute_boundary_neighbors(idx_inner);
      50           0 :     std::vector<coupling::filtering::Quantities<dim>> res(res_offset.size(), 0);
      51             : 
      52           0 :     int tdist = std::min((int)_timeWindowSize, (int)_cycleCounter + 1);
      53             :     // TODO restrict loops to bouncing V_i search window
      54           0 :     for (int offset_t = 0; offset_t > -tdist; offset_t--) {
      55           0 :       unsigned int t = posmod((int)_t + offset_t, _timeWindowSize);
      56           0 :       for (auto idx2 : _innerCellIndices) {
      57           0 :         auto idx2_inner = (tarch::la::Vector<dim, unsigned int>)(idx2.get()) - tarch::la::Vector<dim, unsigned int>{1};
      58           0 :         auto& other = _patchfield(idx2_inner, t);
      59           0 :         double dist = me.distance(other);
      60             : #ifdef NLM_DEBUG
      61             :         std::cout << "dist = " << dist << std::endl;
      62             : #endif
      63           0 :         if (dist < _sigsq)
      64           0 :           cnt_sigsq++;
      65           0 :         if (dist - _sigsq < _hsq)
      66           0 :           cnt_hsq++;
      67           0 :         cnt_total++;
      68           0 :         double weight = std::exp(-std::max(dist - _sigsq, 0.0) / _hsq) * 1 / (-offset_t + 1);
      69             : #ifdef NLM_DEBUG
      70             :         std::cout << "weight = " << weight << std::endl;
      71             : #endif
      72           0 :         for (unsigned int i = 0; i < res_offset.size(); i++) {
      73           0 :           tarch::la::Vector<dim, unsigned int> residx(idx2.get() + res_offset[i]);
      74           0 :           coupling::filtering::Quantities<dim>& q = _flowfield(residx, t);
      75           0 :           res[i] += q * weight;
      76             :         }
      77           0 :         C += weight;
      78             :       }
      79             :     }
      80             : 
      81           0 :     for (unsigned int i = 0; i < res_offset.size(); i++) {
      82           0 :       res[i] = res[i] * (1 / C);
      83           0 :       tarch::la::Vector<dim, int> residx(idx.get() + res_offset[i]);
      84           0 :       int cellidx = coupling::indexing::convertToScalar(CellIndex_T(residx));
      85             : 
      86           0 :       coupling::filtering::JunctorInterface<dim, 2, 1>::_outputCellVectors[0][cellidx]->setMacroscopicMass(res[i][0]);
      87           0 :       tarch::la::Vector<dim, double> momentum;
      88           0 :       for (unsigned int d = 0; d < dim; d++)
      89           0 :         momentum[d] = res[i][d + 1];
      90           0 :       coupling::filtering::JunctorInterface<dim, 2, 1>::_outputCellVectors[0][cellidx]->setMacroscopicMomentum(momentum);
      91             :     }
      92             :   }
      93             :   // dynamic parameter adaption: Goal: cnt/cnt_total should be close to
      94             :   // value_rel
      95           0 :   if (_sigsq_rel >= 0) { // negative value means adaptivity disabled
      96           0 :     double rel = (double)(cnt_sigsq) / cnt_total;
      97           0 :     if (rel < _sigsq_rel * 0.9)
      98           0 :       _sigsq *= 1.05;
      99           0 :     else if (rel > _sigsq_rel * 1.1)
     100           0 :       _sigsq *= 0.95;
     101           0 :     if(_cycleCounter%20 == 0) std::cout << "NLM info: _sigsq is " << _sigsq << std::endl;
     102             :   }
     103           0 :   if (_hsq_rel >= 0) { // negative value means adaptivity disabled
     104           0 :     double rel = (double)(cnt_hsq) / cnt_total;
     105           0 :     if (rel < _hsq_rel * 0.9)
     106           0 :       _hsq *= 1.05;
     107           0 :     else if (rel > _hsq_rel * 1.1)
     108           0 :       _hsq *= 0.95;
     109           0 :     if(_cycleCounter%20 == 0) std::cout << "NLM info: _hsq is " << _hsq << std::endl;
     110             :   }
     111           0 : }
     112             : 
     113             : template <unsigned int dim, coupling::indexing::IndexTrait... scope>
     114             : std::vector<tarch::la::Vector<dim, int>>
     115           0 : coupling::filtering::NLM<dim, scope...>::compute_boundary_neighbors(const tarch::la::Vector<dim, unsigned int>& idx_inner) {
     116           0 :   auto domainSize = CellIndex_T::numberCellsInDomain;
     117           0 :   std::vector<tarch::la::Vector<dim, int>> res_offset;
     118           0 :   res_offset.reserve(1 << dim); // max this many neighbors are possible
     119           0 :   res_offset.push_back(0);      // result for me (i.e. idx) is always needed
     120           0 :   for (unsigned int d = 0; d < dim; d++) {
     121           0 :     tarch::la::Vector<dim, int> offset{0};
     122           0 :     if (idx_inner[d] == 0) {
     123           0 :       offset[d] = -1;
     124           0 :       res_offset.push_back(offset); // Comp res for neighbor cell on boundary surface
     125           0 :     } else if (idx_inner[d] == domainSize[d] - 3) {
     126           0 :       offset[d] = 1;
     127           0 :       res_offset.push_back(offset); // Comp res for neighbor cell on boundary surface
     128             :     }
     129           0 :     if (offset[d] != 0) {
     130           0 :       for (unsigned int d2 = d + 1; d2 < dim; d2++) {
     131           0 :         if (idx_inner[d2] == 0) {
     132           0 :           offset[d2] = -1;
     133           0 :           res_offset.push_back(offset); // Comp res for neighbor cell on boundary edge
     134           0 :         } else if (idx_inner[d2] == domainSize[d2] - 3) {
     135           0 :           offset[d2] = 1;
     136           0 :           res_offset.push_back(offset); // Comp res for neighbor cell on boundary edge
     137             :         }
     138           0 :         if (dim == 3 && d2 == 1 && offset[d2] != 0) {
     139           0 :           if (idx_inner[2] == 0) {
     140           0 :             offset[2] = -1;
     141           0 :             res_offset.push_back(offset); // Comp res for neighbor cell on boundary corner
     142           0 :           } else if (idx_inner[2] == domainSize[2] - 3) {
     143           0 :             offset[2] = 1;
     144           0 :             res_offset.push_back(offset); // Comp res for neighbor cell on boundary corner
     145             :           }
     146           0 :           offset[2] = 0;
     147             :         }
     148           0 :         offset[d2] = 0;
     149             :       }
     150             :     }
     151             :   }
     152           0 :   return res_offset; // (implicit std::move here)
     153           0 : }
     154             : 
     155           0 : template <unsigned int dim, coupling::indexing::IndexTrait... scope> void coupling::filtering::NLM<dim, scope...>::increment_time() {
     156           0 :   _cycleCounter++;
     157           0 :   _t = _cycleCounter % _timeWindowSize;
     158           0 : }