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)
102 0 : std::cout << "NLM info: _sigsq is " << _sigsq << std::endl;
103 : }
104 0 : if (_hsq_rel >= 0) { // negative value means adaptivity disabled
105 0 : double rel = (double)(cnt_hsq) / cnt_total;
106 0 : if (rel < _hsq_rel * 0.9)
107 0 : _hsq *= 1.05;
108 0 : else if (rel > _hsq_rel * 1.1)
109 0 : _hsq *= 0.95;
110 0 : if (_cycleCounter % 20 == 0)
111 0 : std::cout << "NLM info: _hsq is " << _hsq << std::endl;
112 : }
113 0 : }
114 :
115 : template <unsigned int dim, coupling::indexing::IndexTrait... scope>
116 : std::vector<tarch::la::Vector<dim, int>>
117 0 : coupling::filtering::NLM<dim, scope...>::compute_boundary_neighbors(const tarch::la::Vector<dim, unsigned int>& idx_inner) {
118 0 : auto domainSize = CellIndex_T::numberCellsInDomain;
119 0 : std::vector<tarch::la::Vector<dim, int>> res_offset;
120 0 : res_offset.reserve(1 << dim); // max this many neighbors are possible
121 0 : res_offset.push_back(0); // result for me (i.e. idx) is always needed
122 0 : for (unsigned int d = 0; d < dim; d++) {
123 0 : tarch::la::Vector<dim, int> offset{0};
124 0 : if (idx_inner[d] == 0) {
125 0 : offset[d] = -1;
126 0 : res_offset.push_back(offset); // Comp res for neighbor cell on boundary surface
127 0 : } else if (idx_inner[d] == domainSize[d] - 3) {
128 0 : offset[d] = 1;
129 0 : res_offset.push_back(offset); // Comp res for neighbor cell on boundary surface
130 : }
131 0 : if (offset[d] != 0) {
132 0 : for (unsigned int d2 = d + 1; d2 < dim; d2++) {
133 0 : if (idx_inner[d2] == 0) {
134 0 : offset[d2] = -1;
135 0 : res_offset.push_back(offset); // Comp res for neighbor cell on boundary edge
136 0 : } else if (idx_inner[d2] == domainSize[d2] - 3) {
137 0 : offset[d2] = 1;
138 0 : res_offset.push_back(offset); // Comp res for neighbor cell on boundary edge
139 : }
140 0 : if (dim == 3 && d2 == 1 && offset[d2] != 0) {
141 0 : if (idx_inner[2] == 0) {
142 0 : offset[2] = -1;
143 0 : res_offset.push_back(offset); // Comp res for neighbor cell on boundary corner
144 0 : } else if (idx_inner[2] == domainSize[2] - 3) {
145 0 : offset[2] = 1;
146 0 : res_offset.push_back(offset); // Comp res for neighbor cell on boundary corner
147 : }
148 0 : offset[2] = 0;
149 : }
150 0 : offset[d2] = 0;
151 : }
152 : }
153 : }
154 0 : return res_offset; // (implicit std::move here)
155 0 : }
156 :
157 0 : template <unsigned int dim, coupling::indexing::IndexTrait... scope> void coupling::filtering::NLM<dim, scope...>::increment_time() {
158 0 : _cycleCounter++;
159 0 : _t = _cycleCounter % _timeWindowSize;
160 0 : }
|
