Line data Source code
1 : // Copyright (C) 2015 Technische Universitaet Muenchen
2 : // This file is part of the Mamico project. For conditions of distribution
3 : // and use, please see the copyright notice in Mamico's main folder, or at
4 : // www5.in.tum.de/mamico
5 :
6 : template <class LinkedCell, unsigned int dim>
7 0 : coupling::transferstrategies::TransferStrategy4NieCoupling<LinkedCell, dim>::TransferStrategy4NieCoupling(
8 : coupling::interface::MDSolverInterface<LinkedCell, dim>* const mdSolverInterface, unsigned int numberMDSteps, double shiftTimestep,
9 : tarch::la::Vector<2 * dim, bool> massFluxBoundary, bool energyCorrection, double density)
10 0 : : coupling::transferstrategies::TransferStrategy<LinkedCell, dim>(mdSolverInterface), _massMapping(mdSolverInterface), _momentumMapping(mdSolverInterface),
11 0 : _oldSolution(NULL), _newSolution(NULL), _numberMDSteps(numberMDSteps), _shiftTimestep(shiftTimestep), _numberLocalCells(I02::linearNumberCellsInDomain),
12 0 : _excessMass(NULL), _massFluxBoundary(massFluxBoundary), _enableEnergyCorrection(energyCorrection), _density(density) {
13 : // allocate arrays for old and new solution (required for linear time
14 : // interpolation) and init. them with zero vectors
15 0 : _oldSolution = new tarch::la::Vector<dim, double>[_numberLocalCells];
16 0 : _newSolution = new tarch::la::Vector<dim, double>[_numberLocalCells];
17 0 : _excessMass = new double[_numberLocalCells];
18 0 : if ((_oldSolution == NULL) || (_newSolution == NULL) || (_excessMass == NULL)) {
19 0 : std::cout << "ERROR "
20 : "coupling::transferstrategies::TransferStrategy4NieCoupling::"
21 : "TransferStrategy4NieCoupling(...): _ptr==NULL!"
22 0 : << std::endl;
23 0 : exit(EXIT_FAILURE);
24 : }
25 0 : for (unsigned int i = 0; i < _numberLocalCells; i++) {
26 0 : _oldSolution[i] = tarch::la::Vector<dim, double>(0.0);
27 0 : _newSolution[i] = tarch::la::Vector<dim, double>(0.0);
28 0 : _excessMass[i] = 0.0;
29 : }
30 0 : }
31 :
32 0 : template <class LinkedCell, unsigned int dim> coupling::transferstrategies::TransferStrategy4NieCoupling<LinkedCell, dim>::~TransferStrategy4NieCoupling() {
33 : // delete arrays
34 0 : if (_oldSolution != NULL) {
35 0 : delete[] _oldSolution;
36 0 : _oldSolution = NULL;
37 : }
38 0 : if (_newSolution != NULL) {
39 0 : delete[] _newSolution;
40 0 : _newSolution = NULL;
41 : }
42 0 : if (_excessMass != NULL) {
43 0 : delete[] _excessMass;
44 0 : _excessMass = NULL;
45 : }
46 0 : }
47 :
48 : template <class LinkedCell, unsigned int dim>
49 0 : void coupling::transferstrategies::TransferStrategy4NieCoupling<LinkedCell, dim>::beginProcessInnerCouplingCellsBeforeReceivingMacroscopicSolverData() {
50 : // reset time counter (always counts from 0 to _numberMDSteps)
51 0 : _timestepCounter = 0;
52 : //_totalMass=0;
53 : //_cellCount=0;
54 : // copy new solution to old solution
55 0 : for (unsigned int i = 0; i < _numberLocalCells; i++) {
56 0 : _oldSolution[i] = _newSolution[i];
57 : }
58 0 : }
59 :
60 : template <class LinkedCell, unsigned int dim>
61 0 : void coupling::transferstrategies::TransferStrategy4NieCoupling<LinkedCell, dim>::processInnerCouplingCellBeforeReceivingMacroscopicSolverData(
62 : coupling::datastructures::CouplingCellWithLinkedCells<LinkedCell, dim>& cell, I02 index) {
63 : // backup old mass that could not be inserted into cells
64 0 : _excessMass[index.get()] = cell.getMicroscopicMass();
65 :
66 : // TODO: remove temporary debug output
67 : // if(_excessMass[index.get()] > 0) std::cout << "TransferStrategy4NieCoupling:
68 : // _excessMass[" << index << "]: " << _excessMass[index.get()] << std::endl;
69 : // TODO: make threshold customizable
70 : // const unsigned int threshold = 3;
71 : /*if(_excessMass[index.get()] > threshold) {
72 : std::cout << "TransferStrategy4NieCoupling: ERROR: Excess mass above
73 : threshold (= " << threshold << "): _excessMass[" << index << "] = " <<
74 : _excessMass[index.get()] << std::endl; exit(EXIT_FAILURE);
75 : }*/
76 :
77 : // reset cell-mass (would be overwritten anyway; but just to be on safe side)
78 0 : cell.setMicroscopicMass(0.0);
79 0 : }
80 :
81 : template <class LinkedCell, unsigned int dim>
82 0 : void coupling::transferstrategies::TransferStrategy4NieCoupling<LinkedCell, dim>::processInnerCouplingCellAfterReceivingMacroscopicSolverData(
83 : coupling::datastructures::CouplingCellWithLinkedCells<LinkedCell, dim>& cell, I02 index) {
84 :
85 : // convert momentum to velocity values for cont->MD transfer
86 0 : if (cell.getMicroscopicMass() == 0.0) {
87 0 : cell.setMicroscopicMomentum(tarch::la::Vector<dim, double>(0.0));
88 : } else {
89 0 : cell.setMicroscopicMomentum((1.0 / cell.getMicroscopicMass()) * cell.getMicroscopicMomentum());
90 : }
91 :
92 : // write velocity value to new solution array and interpolate first solution
93 : // (= value of old solution)
94 0 : _newSolution[index.get()] = cell.getMicroscopicMomentum();
95 0 : cell.setMicroscopicMomentum((1.0 - _shiftTimestep) * _oldSolution[index.get()] + _shiftTimestep * _newSolution[index.get()]);
96 :
97 0 : const tarch::la::Vector<dim, double> velocity = (0.5 - _shiftTimestep) * _oldSolution[index.get()] + (_shiftTimestep + 0.5) * _newSolution[index.get()];
98 :
99 : // for mass transfer: compute mass flux if required and add excess mass
100 0 : const double massFlux = computeMassFlux(cell.getMicroscopicMass(), velocity, index) + _excessMass[index.get()];
101 0 : cell.setMicroscopicMass(massFlux);
102 :
103 : // reset macroscopic solver values (from averaging)
104 0 : cell.setMacroscopicMass(0.0);
105 0 : cell.setMacroscopicMomentum(tarch::la::Vector<dim, double>(0.0));
106 0 : }
107 :
108 : // template<class LinkedCell,unsigned int dim>
109 : // void
110 : // coupling::transferstrategies::TransferStrategy4NieCoupling<LinkedCell,dim>::
111 : // beginProcessInnerCouplingCellsBeforeSendingMDSolverData(){
112 : // }
113 :
114 : template <class LinkedCell, unsigned int dim>
115 0 : void coupling::transferstrategies::TransferStrategy4NieCoupling<LinkedCell, dim>::processInnerCouplingCellBeforeSendingMDSolverData(
116 : coupling::datastructures::CouplingCellWithLinkedCells<LinkedCell, dim>& cell, I02 index) {
117 : #if (COUPLING_MD_DEBUG == COUPLING_MD_YES)
118 : std::cout << "processInnerCouplingCellBeforeSendingMDSolverData(): Data "
119 : "before averaging from cell "
120 : << index << "Mass: " << cell.getMacroscopicMass() << " , momentum: " << cell.getMacroscopicMomentum() << std::endl;
121 : #endif
122 :
123 : // average sampled data from MD
124 0 : if (_timestepCounter != 0) {
125 0 : const double mass = cell.getMacroscopicMass() / ((double)_timestepCounter);
126 0 : double momentumFactor = (1.0 / ((double)_timestepCounter));
127 :
128 0 : if (_enableEnergyCorrection) {
129 0 : const double mass_macro = IDXS.getCouplingCellSize()[0] * IDXS.getCouplingCellSize()[1] * IDXS.getCouplingCellSize()[2] * _density;
130 : // computes energy conservation correction factor
131 : // based on assumption mass/2*vel_micro^2 = mass_macro/2*vel_macro^2
132 : // Note: In MD-60 scenario, energy_correction factor is on average ca. 0.9998
133 0 : const double energy_correction = sqrt(mass / mass_macro);
134 0 : momentumFactor *= energy_correction;
135 : }
136 :
137 0 : const tarch::la::Vector<dim, double> momentum = momentumFactor * cell.getMacroscopicMomentum();
138 0 : cell.setMacroscopicMass(mass);
139 0 : cell.setMacroscopicMomentum(momentum);
140 : } else {
141 0 : cell.setMacroscopicMass(0.0);
142 0 : cell.setMacroscopicMomentum(tarch::la::Vector<dim, double>(0.0));
143 : }
144 : #if (COUPLING_MD_DEBUG == COUPLING_MD_YES)
145 : std::cout << "Time step counter=" << _timestepCounter << std::endl;
146 : std::cout << "processInnerCouplingCellBeforeSendingMDSolverData(): Data "
147 : "from cell "
148 : << index << "Mass: " << cell.getMacroscopicMass() << " , momentum: " << cell.getMacroscopicMomentum() << std::endl;
149 : #endif
150 :
151 : // cell.iterateConstCells(_massMapping);
152 : // _totalMass += _massMapping.getMass();
153 : // _totalMass += cell.getMicroscopicMass();
154 0 : }
155 :
156 : // template<class LinkedCell,unsigned int dim>
157 : // void
158 : // coupling::transferstrategies::TransferStrategy4NieCoupling<LinkedCell,dim>::
159 : // endProcessInnerCouplingCellsBeforeSendingMDSolverData(){
160 : // std::cout << "Avg massFluxBoundary MD density: " <<
161 : // _totalMass/(_cellCount*15.625) << std::endl;
162 : //}
163 :
164 : template <class LinkedCell, unsigned int dim>
165 0 : void coupling::transferstrategies::TransferStrategy4NieCoupling<LinkedCell, dim>::beginProcessInnerCouplingCellsAfterMDTimestep() {
166 : // increment time step counter
167 0 : _timestepCounter++;
168 0 : }
169 :
170 : template <class LinkedCell, unsigned int dim>
171 0 : void coupling::transferstrategies::TransferStrategy4NieCoupling<LinkedCell, dim>::processInnerCouplingCellAfterMDTimestep(
172 : coupling::datastructures::CouplingCellWithLinkedCells<LinkedCell, dim>& cell, I02 index) {
173 : // sampling: add up mass and momentum in macroscopic buffers
174 0 : cell.iterateConstCells(_massMapping);
175 0 : const double mass = _massMapping.getMass();
176 : //_totalMass += mass;
177 : //_cellCount++;
178 0 : cell.iterateConstCells(_momentumMapping);
179 0 : const tarch::la::Vector<dim, double> momentum(_momentumMapping.getMomentum());
180 0 : cell.addMacroscopicMass(mass);
181 0 : cell.addMacroscopicMomentum(momentum);
182 :
183 : // interpolation of velocity value from old and new solution
184 : // we shift the factor by _shiftTimestep time steps and inter- or extrapolate
185 : // velocities. A CFD solver will deliver values at t and t+dt and we want to
186 : // advance the simulation potentially differently
187 0 : const double factor = _shiftTimestep + ((double)_timestepCounter) / ((double)_numberMDSteps);
188 : #if (COUPLING_MD_DEBUG == COUPLING_MD_YES)
189 : if (factor < _shiftTimestep || factor > 1.0 + _shiftTimestep) {
190 : std::cout << "ERROR "
191 : "coupling::transferstrategies::TransferStrategy4NieCoupling<"
192 : "LinkedCell,dim>::processInnerCouplingCellAfterMDTimestep(."
193 : ".): factor out of range!"
194 : << std::endl;
195 : exit(EXIT_FAILURE);
196 : }
197 : #endif
198 0 : const tarch::la::Vector<dim, double> velocity = (1.0 - factor) * _oldSolution[index.get()] + factor * _newSolution[index.get()];
199 0 : cell.setMicroscopicMomentum(velocity);
200 0 : }
201 :
202 : template <class LinkedCell, unsigned int dim>
203 0 : double coupling::transferstrategies::TransferStrategy4NieCoupling<LinkedCell, dim>::computeMassFlux(const double& mass,
204 : const tarch::la::Vector<dim, double>& velocity,
205 : const I01 index) {
206 : // meshsize and volume of cell
207 0 : const tarch::la::Vector<dim, double> meshsize(IDXS.getCouplingCellSize());
208 : double volume = 1.0;
209 0 : for (unsigned int d = 0; d < dim; d++) {
210 0 : volume = volume * meshsize[d];
211 : }
212 :
213 : // mass density of cell
214 0 : const double density = mass / volume;
215 : // coupling time interval
216 0 : const double couplingTimeInterval = coupling::transferstrategies::TransferStrategy<LinkedCell, dim>::_mdSolverInterface->getDt() * _numberMDSteps;
217 :
218 0 : double massFlux = 0.0;
219 :
220 : // compute normal vector for this cell, checking whether there is mass flux or
221 : // not
222 0 : tarch::la::Vector<dim, int> n(0);
223 0 : for (unsigned int d = 0; d < dim; d++) {
224 : // compute normal vector (pointing inwards MD domain)
225 0 : if ((index.get()[d] == 1) && _massFluxBoundary[2 * d]) {
226 0 : n[d] = 1; /*_totalMass+=mass;_cellCount++;*/
227 0 : } else if ((index.get()[d] == (int)I09::numberCellsInDomain[d]) && _massFluxBoundary[2 * d + 1]) {
228 0 : n[d] = -1; /*_totalMass+=mass;_cellCount++;*/
229 : }
230 : }
231 :
232 : // compute mass flux for each boundary
233 0 : for (unsigned int d = 0; d < dim; d++) {
234 : double surface = 1.0;
235 0 : for (unsigned int e = 0; e < d; e++) {
236 0 : surface = surface * meshsize[e];
237 : }
238 0 : for (unsigned int e = d + 1; e < dim; e++) {
239 0 : surface = surface * meshsize[e];
240 : }
241 :
242 0 : massFlux += n[d] * surface * density * velocity[d] * couplingTimeInterval;
243 : }
244 0 : return massFlux;
245 : }
|
