Line data Source code
1 : // EnergyPlus, Copyright (c) 1996-2025, The Board of Trustees of the University of Illinois,
2 : // The Regents of the University of California, through Lawrence Berkeley National Laboratory
3 : // (subject to receipt of any required approvals from the U.S. Dept. of Energy), Oak Ridge
4 : // National Laboratory, managed by UT-Battelle, Alliance for Sustainable Energy, LLC, and other
5 : // contributors. All rights reserved.
6 : //
7 : // NOTICE: This Software was developed under funding from the U.S. Department of Energy and the
8 : // U.S. Government consequently retains certain rights. As such, the U.S. Government has been
9 : // granted for itself and others acting on its behalf a paid-up, nonexclusive, irrevocable,
10 : // worldwide license in the Software to reproduce, distribute copies to the public, prepare
11 : // derivative works, and perform publicly and display publicly, and to permit others to do so.
12 : //
13 : // Redistribution and use in source and binary forms, with or without modification, are permitted
14 : // provided that the following conditions are met:
15 : //
16 : // (1) Redistributions of source code must retain the above copyright notice, this list of
17 : // conditions and the following disclaimer.
18 : //
19 : // (2) Redistributions in binary form must reproduce the above copyright notice, this list of
20 : // conditions and the following disclaimer in the documentation and/or other materials
21 : // provided with the distribution.
22 : //
23 : // (3) Neither the name of the University of California, Lawrence Berkeley National Laboratory,
24 : // the University of Illinois, U.S. Dept. of Energy nor the names of its contributors may be
25 : // used to endorse or promote products derived from this software without specific prior
26 : // written permission.
27 : //
28 : // (4) Use of EnergyPlus(TM) Name. If Licensee (i) distributes the software in stand-alone form
29 : // without changes from the version obtained under this License, or (ii) Licensee makes a
30 : // reference solely to the software portion of its product, Licensee must refer to the
31 : // software as "EnergyPlus version X" software, where "X" is the version number Licensee
32 : // obtained under this License and may not use a different name for the software. Except as
33 : // specifically required in this Section (4), Licensee shall not use in a company name, a
34 : // product name, in advertising, publicity, or other promotional activities any name, trade
35 : // name, trademark, logo, or other designation of "EnergyPlus", "E+", "e+" or confusingly
36 : // similar designation, without the U.S. Department of Energy's prior written consent.
37 : //
38 : // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
39 : // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
40 : // AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
41 : // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
42 : // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
43 : // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
44 : // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
45 : // OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
46 : // POSSIBILITY OF SUCH DAMAGE.
47 :
48 : // C++ Headers
49 : #include <cassert>
50 : #include <cmath>
51 :
52 : // ObjexxFCL Headers
53 : #include <ObjexxFCL/Array.functions.hh>
54 : #include <ObjexxFCL/ArrayS.functions.hh>
55 : #include <ObjexxFCL/Fmath.hh>
56 :
57 : // EnergyPlus Headers
58 : #include <EnergyPlus/Construction.hh>
59 : #include <EnergyPlus/Data/EnergyPlusData.hh>
60 : #include <EnergyPlus/DataHeatBalSurface.hh>
61 : #include <EnergyPlus/DataHeatBalance.hh>
62 : #include <EnergyPlus/DataIPShortCuts.hh>
63 : #include <EnergyPlus/DataSurfaces.hh>
64 : #include <EnergyPlus/DataSystemVariables.hh>
65 : #include <EnergyPlus/DataViewFactorInformation.hh>
66 : #include <EnergyPlus/DisplayRoutines.hh>
67 : #include <EnergyPlus/General.hh>
68 : #include <EnergyPlus/HeatBalanceIntRadExchange.hh>
69 : #include <EnergyPlus/InputProcessing/InputProcessor.hh>
70 : #include <EnergyPlus/Material.hh>
71 : #include <EnergyPlus/UtilityRoutines.hh>
72 : #include <EnergyPlus/WindowEquivalentLayer.hh>
73 :
74 : namespace EnergyPlus {
75 :
76 : namespace HeatBalanceIntRadExchange {
77 : // Module containing the routines dealing with the interior radiant exchange
78 : // between surfaces.
79 :
80 : // MODULE INFORMATION:
81 : // AUTHOR Rick Strand
82 : // DATE WRITTEN September 2000
83 : // MODIFIED Aug 2001, FW: recalculate ScriptF for a zone if window interior
84 : // shade/blind status is different from previous time step. This is
85 : // because ScriptF, which is used to calculate interior LW
86 : // exchange between surfaces, depends on inside surface emissivities,
87 : // which, for a window, depends on whether or not an interior
88 : // shade or blind is in place.
89 :
90 : // PURPOSE OF THIS MODULE:
91 : // Part of the heat balance modularization/re-engineering. Purpose of this
92 : // module is to replace the MRT with RBAL method of modeling radiant exchange
93 : // between interior surfaces.
94 :
95 : // METHODOLOGY EMPLOYED:
96 : // Standard EnergyPlus methodology
97 :
98 : // REFERENCES:
99 : // ASHRAE Loads Toolkit "Script F" routines by Curt Pedersen
100 : // Hottel, H.C., and A.F. Sarofim. "Radiative Transfer" (mainly chapter 3),
101 : // McGraw-Hill, Inc., New York, 1967.
102 :
103 1377810 : void CalcInteriorRadExchange(EnergyPlusData &state,
104 : Array1S<Real64> const SurfaceTemp, // Current surface temperatures
105 : int const SurfIterations, // Number of iterations in calling subroutine
106 : Array1D<Real64> &NetLWRadToSurf, // Net long wavelength radiant exchange from other surfaces
107 : ObjexxFCL::Optional_int_const ZoneToResimulate, // if passed in, then only calculate for this zone
108 : #ifdef EP_Count_Calls
109 : std::string_view const CalledFrom)
110 : #else
111 : [[maybe_unused]] std::string_view const CalledFrom)
112 : #endif
113 : {
114 :
115 : // SUBROUTINE INFORMATION:
116 : // AUTHOR Rick Strand
117 : // DATE WRITTEN September 2000
118 : // MODIFIED 6/18/01, FCW: calculate IR on windows
119 : // Jan 2002, FCW: add blinds with movable slats
120 : // Sep 2011 LKL/BG - resimulate only zones needing it for Radiant systems
121 :
122 : // PURPOSE OF THIS SUBROUTINE:
123 : // Determines the interior radiant exchange between surfaces using
124 : // Hottel's ScriptF method for the grey interchange between surfaces
125 : // in an enclosure.
126 :
127 : // REFERENCES:
128 : // Hottel, H. C. and A. F. Sarofim, Radiative Transfer, Ch 3, McGraw Hill, 1967.
129 :
130 : // Types
131 : typedef Array1D<Real64>::size_type size_type;
132 :
133 : bool IntShadeOrBlindStatusChanged; // True if status of interior shade or blind on at least
134 : // one window in a zone has changed from previous time step
135 :
136 1377810 : auto &SurfaceTempRad = state.dataHeatBalIntRadExchg->SurfaceTempRad;
137 1377810 : auto &SurfaceTempInKto4th = state.dataHeatBalIntRadExchg->SurfaceTempInKto4th;
138 1377810 : auto &SurfaceEmiss = state.dataHeatBalIntRadExchg->SurfaceEmiss;
139 :
140 1377810 : if (state.dataHeatBalIntRadExchg->CalcInteriorRadExchangefirstTime) {
141 114 : SurfaceTempRad.allocate(state.dataHeatBalIntRadExchg->MaxNumOfRadEnclosureSurfs);
142 114 : SurfaceTempInKto4th.allocate(state.dataHeatBalIntRadExchg->MaxNumOfRadEnclosureSurfs);
143 114 : SurfaceEmiss.allocate(state.dataHeatBalIntRadExchg->MaxNumOfRadEnclosureSurfs);
144 114 : state.dataHeatBalIntRadExchg->CalcInteriorRadExchangefirstTime = false;
145 114 : if (state.dataSysVars->DeveloperFlag) {
146 0 : DisplayString(state, " OMP turned off, HBIRE loop executed in serial");
147 : }
148 : }
149 :
150 1377810 : if (state.dataGlobal->KickOffSimulation || state.dataGlobal->KickOffSizing) return;
151 :
152 1362285 : bool const PartialResimulate(present(ZoneToResimulate));
153 :
154 : #ifdef EP_Count_Calls
155 : if (!PartialResimulate) {
156 : ++state.dataTimingsData->NumIntRadExchange_Calls;
157 : } else {
158 : ++state.dataTimingsData->NumIntRadExchangeZ_Calls;
159 : }
160 : if (CalledFrom.empty()) {
161 : // do nothing
162 : } else if (CalledFrom == "Main") {
163 : ++state.dataTimingsData->NumIntRadExchangeMain_Calls;
164 : } else if (CalledFrom == "Outside") {
165 : ++state.dataTimingsData->NumIntRadExchangeOSurf_Calls;
166 : } else if (CalledFrom == "Inside") {
167 : ++state.dataTimingsData->NumIntRadExchangeISurf_Calls;
168 : }
169 : #endif
170 :
171 1362285 : int startEnclosure = 1;
172 1362285 : int endEnclosure = state.dataViewFactor->NumOfRadiantEnclosures;
173 1362285 : if (PartialResimulate) {
174 : // ToDo: For now, use min and max enclosure numbers associated with this zone, this could include unrelated enclosures
175 0 : startEnclosure = state.dataHeatBal->Zone(ZoneToResimulate).zoneRadEnclosureFirst;
176 0 : endEnclosure = state.dataHeatBal->Zone(ZoneToResimulate).zoneRadEnclosureLast;
177 0 : for (int enclosureNum = startEnclosure; enclosureNum <= endEnclosure; ++enclosureNum) {
178 0 : auto const &enclosure = state.dataViewFactor->EnclRadInfo(enclosureNum);
179 0 : for (int i : enclosure.SurfacePtr) {
180 0 : NetLWRadToSurf(i) = 0.0;
181 0 : state.dataSurface->SurfWinIRfromParentZone(i) = 0.0;
182 : }
183 : }
184 : } else {
185 1362285 : NetLWRadToSurf = 0.0;
186 13388294 : for (int SurfNum = 1; SurfNum <= state.dataSurface->TotSurfaces; SurfNum++)
187 12026009 : state.dataSurface->SurfWinIRfromParentZone(SurfNum) = 0.0;
188 : }
189 :
190 3354607 : for (int enclosureNum = startEnclosure; enclosureNum <= endEnclosure; ++enclosureNum) {
191 :
192 1992322 : auto &zone_info = state.dataViewFactor->EnclRadInfo(enclosureNum);
193 1992322 : auto &zone_ScriptF = zone_info.ScriptF; // Tuned Transposed
194 1992322 : int const n_zone_Surfaces = zone_info.NumOfSurfaces;
195 1992322 : size_type const s_zone_Surfaces = n_zone_Surfaces;
196 :
197 : // Calculate ScriptF if first time step in environment and surface heat-balance iterations not yet started;
198 : // recalculate ScriptF if status of window interior shades or blinds has changed from
199 : // previous time step. This recalculation is required since ScriptF depends on the inside
200 : // emissivity of the inside surfaces, which, for windows, is (1) the emissivity of the
201 : // inside face of the inside glass layer if there is no interior shade/blind, or (2) the effective
202 : // emissivity of the shade/blind if the shade/blind is in place. (The "effective emissivity"
203 : // in this case is (1) the shade/blind emissivity if the shade/blind IR transmittance is zero,
204 : // or (2) a weighted average of the shade/blind emissivity and inside glass emissivity if the
205 : // shade/blind IR transmittance is not zero (which is sometimes the case for a "shade" and
206 : // usually the case for a blind). It assumed for switchable glazing that the inside surface
207 : // emissivity does not change if the glazing is switched on or off.
208 :
209 : // Determine if status of interior shade/blind on one or more windows in the zone has changed
210 : // from previous time step. Also make a check for any changes in interior movable insulation.
211 :
212 1992322 : if (SurfIterations == 0) {
213 :
214 : bool IntMovInsulChanged; // True if the status of interior movable insulation has changed
215 :
216 1063282 : IntShadeOrBlindStatusChanged = false;
217 1063282 : IntMovInsulChanged = false;
218 :
219 1063282 : if (!state.dataGlobal->BeginEnvrnFlag) { // Check for change in shade/blind status
220 7259987 : for (int const SurfNum : zone_info.SurfacePtr) {
221 6197755 : if (IntShadeOrBlindStatusChanged || IntMovInsulChanged)
222 : break; // Need only check if one window's status or one movable insulation status has changed
223 6197755 : if (state.dataConstruction->Construct(state.dataSurface->Surface(SurfNum).Construction).TypeIsWindow) {
224 264364 : DataSurfaces::WinShadingType ShadeFlag = state.dataSurface->SurfWinShadingFlag(SurfNum);
225 264364 : DataSurfaces::WinShadingType ShadeFlagPrev = state.dataSurface->SurfWinExtIntShadePrevTS(SurfNum);
226 264364 : if (ShadeFlagPrev != ShadeFlag && (ANY_INTERIOR_SHADE_BLIND(ShadeFlagPrev) || ANY_INTERIOR_SHADE_BLIND(ShadeFlag)))
227 0 : IntShadeOrBlindStatusChanged = true;
228 271407 : if (state.dataSurface->SurfWinWindowModelType(SurfNum) == DataSurfaces::WindowModel::EQL &&
229 : state.dataWindowEquivLayer
230 7043 : ->CFS(state.dataConstruction->Construct(state.dataSurface->Surface(SurfNum).Construction).EQLConsPtr)
231 7043 : .ISControlled) {
232 2874 : IntShadeOrBlindStatusChanged = true;
233 : }
234 : } else {
235 5933391 : if (state.dataSurface->AnyMovableInsulation) UpdateMovableInsulationFlag(state, IntMovInsulChanged, SurfNum);
236 : }
237 : }
238 : }
239 :
240 2123690 : if (IntShadeOrBlindStatusChanged || IntMovInsulChanged ||
241 1060408 : state.dataGlobal->BeginEnvrnFlag) { // Calc inside surface emissivities for this time step
242 30138 : for (int ZoneSurfNum = 1; ZoneSurfNum <= n_zone_Surfaces; ++ZoneSurfNum) {
243 26214 : int const SurfNum = zone_info.SurfacePtr(ZoneSurfNum);
244 26214 : int const ConstrNum = state.dataSurface->Surface(SurfNum).Construction;
245 26214 : zone_info.Emissivity(ZoneSurfNum) = state.dataHeatBalSurf->SurfAbsThermalInt(SurfNum);
246 29445 : if (state.dataConstruction->Construct(ConstrNum).TypeIsWindow &&
247 3231 : ANY_INTERIOR_SHADE_BLIND(state.dataSurface->SurfWinShadingFlag(SurfNum))) {
248 0 : zone_info.Emissivity(ZoneSurfNum) = state.dataHeatBalSurf->SurfAbsThermalInt(SurfNum);
249 : }
250 29103 : if (state.dataSurface->SurfWinWindowModelType(SurfNum) == DataSurfaces::WindowModel::EQL &&
251 2889 : state.dataWindowEquivLayer->CFS(state.dataConstruction->Construct(ConstrNum).EQLConsPtr).ISControlled) {
252 2880 : zone_info.Emissivity(ZoneSurfNum) = WindowEquivalentLayer::EQLWindowInsideEffectiveEmiss(state, ConstrNum);
253 : }
254 : }
255 :
256 3924 : if (state.dataHeatBalIntRadExchg->CarrollMethod) {
257 0 : CalcFp(n_zone_Surfaces, zone_info.Emissivity, zone_info.FMRT, zone_info.Fp);
258 : } else {
259 3924 : CalcScriptF(state, n_zone_Surfaces, zone_info.Area, zone_info.F, zone_info.Emissivity, zone_ScriptF);
260 : // precalc - multiply by StefanBoltzmannConstant
261 3924 : zone_ScriptF *= Constant::StefanBoltzmann;
262 : }
263 : }
264 :
265 : } // End of check if SurfIterations = 0
266 :
267 : // Set surface emissivities and temperatures
268 : // Also, for Carroll method, calculate numerators and denominators of radiant temperature
269 1992322 : Real64 CarrollMRTNumerator(0.0);
270 1992322 : Real64 CarrollMRTDenominator(0.0);
271 : Real64 CarrollMRTInKTo4th; // Carroll MRT
272 13788155 : for (size_type ZoneSurfNum = 0; ZoneSurfNum < s_zone_Surfaces; ++ZoneSurfNum) {
273 11795833 : int const SurfNum = zone_info.SurfacePtr[ZoneSurfNum];
274 11795833 : auto const &surf = state.dataSurface->Surface(SurfNum);
275 11795833 : auto const &surfWindow = state.dataSurface->SurfaceWindow(SurfNum);
276 11795833 : int const constrNum = surf.Construction;
277 11795833 : auto const &construct = state.dataConstruction->Construct(constrNum);
278 11795833 : if (construct.WindowTypeEQL) {
279 25151 : SurfaceTempRad[ZoneSurfNum] = state.dataSurface->SurfWinEffInsSurfTemp(SurfNum);
280 25151 : SurfaceEmiss[ZoneSurfNum] = WindowEquivalentLayer::EQLWindowInsideEffectiveEmiss(state, constrNum);
281 11770682 : } else if (construct.WindowTypeBSDF && state.dataSurface->SurfWinShadingFlag(SurfNum) == DataSurfaces::WinShadingType::IntShade) {
282 0 : auto &surfShade = state.dataSurface->surfShades(SurfNum);
283 0 : SurfaceTempRad[ZoneSurfNum] = state.dataSurface->SurfWinEffInsSurfTemp(SurfNum);
284 0 : SurfaceEmiss[ZoneSurfNum] = surfShade.effShadeEmi + surfShade.effGlassEmi;
285 11770682 : } else if (construct.WindowTypeBSDF) {
286 10 : SurfaceTempRad[ZoneSurfNum] = state.dataSurface->SurfWinEffInsSurfTemp(SurfNum);
287 10 : SurfaceEmiss[ZoneSurfNum] = construct.InsideAbsorpThermal;
288 11770672 : } else if (construct.TypeIsWindow && surf.OriginalClass != DataSurfaces::SurfaceClass::TDD_Diffuser) {
289 1048120 : if (SurfIterations == 0 && NOT_SHADED(state.dataSurface->SurfWinShadingFlag(SurfNum))) {
290 : // If the window is bare this TS and it is the first time through we use the previous TS glass
291 : // temperature whether or not the window was shaded in the previous TS. If the window was shaded
292 : // the previous time step this temperature is a better starting value than the shade temperature.
293 257660 : SurfaceTempRad[ZoneSurfNum] = surfWindow.thetaFace[2 * construct.TotGlassLayers] - Constant::Kelvin;
294 257660 : SurfaceEmiss[ZoneSurfNum] = construct.InsideAbsorpThermal;
295 : // For windows with an interior shade or blind an effective inside surface temp
296 : // and emiss is used here that is a weighted combination of shade/blind and glass temp and emiss.
297 266400 : } else if (ANY_INTERIOR_SHADE_BLIND(state.dataSurface->SurfWinShadingFlag(SurfNum))) {
298 0 : SurfaceTempRad[ZoneSurfNum] = state.dataSurface->SurfWinEffInsSurfTemp(SurfNum);
299 0 : SurfaceEmiss[ZoneSurfNum] = state.dataHeatBalSurf->SurfAbsThermalInt(SurfNum);
300 : } else {
301 266400 : SurfaceTempRad[ZoneSurfNum] = SurfaceTemp(SurfNum);
302 266400 : SurfaceEmiss[ZoneSurfNum] = construct.InsideAbsorpThermal;
303 : }
304 : } else {
305 11246612 : SurfaceTempRad[ZoneSurfNum] = SurfaceTemp(SurfNum);
306 11246612 : SurfaceEmiss[ZoneSurfNum] = construct.InsideAbsorpThermal;
307 : }
308 11795833 : if (state.dataHeatBalIntRadExchg->CarrollMethod) {
309 0 : CarrollMRTNumerator += SurfaceTempRad[ZoneSurfNum] * zone_info.Fp[ZoneSurfNum] * zone_info.Area[ZoneSurfNum];
310 0 : CarrollMRTDenominator += zone_info.Fp[ZoneSurfNum] * zone_info.Area[ZoneSurfNum];
311 : }
312 11795833 : SurfaceTempInKto4th[ZoneSurfNum] = pow_4(SurfaceTempRad[ZoneSurfNum] + Constant::Kelvin);
313 : }
314 :
315 1992322 : if (state.dataHeatBalIntRadExchg->CarrollMethod) {
316 0 : if (CarrollMRTDenominator > 0.0) {
317 0 : CarrollMRTInKTo4th = pow_4(CarrollMRTNumerator / CarrollMRTDenominator + Constant::Kelvin);
318 : } else {
319 : // Likely only one surface in this enclosure
320 0 : CarrollMRTInKTo4th = 293.15; // arbitrary value, IR will be zero
321 : }
322 : // These are the money loops
323 0 : for (size_type RecZoneSurfNum = 0; RecZoneSurfNum < s_zone_Surfaces; ++RecZoneSurfNum) {
324 0 : int const RecSurfNum = zone_info.SurfacePtr[RecZoneSurfNum];
325 0 : int const ConstrNumRec = state.dataSurface->Surface(RecSurfNum).Construction;
326 0 : auto const &rec_construct = state.dataConstruction->Construct(ConstrNumRec);
327 0 : auto &netLWRadToRecSurf = NetLWRadToSurf(RecSurfNum);
328 0 : if (rec_construct.TypeIsWindow) {
329 : // auto& rec_surface_window(SurfaceWindow(RecSurfNum));
330 0 : Real64 CarrollMRTInKTo4thWin = CarrollMRTInKTo4th; // arbitrary value, IR will be zero
331 0 : Real64 CarrollMRTNumeratorWin(0.0);
332 0 : Real64 CarrollMRTDenominatorWin(0.0);
333 0 : for (size_type SendZoneSurfNum = 0; SendZoneSurfNum < s_zone_Surfaces; ++SendZoneSurfNum) {
334 0 : if (SendZoneSurfNum != RecZoneSurfNum) {
335 0 : CarrollMRTNumeratorWin +=
336 0 : SurfaceTempRad[SendZoneSurfNum] * zone_info.Fp[SendZoneSurfNum] * zone_info.Area[SendZoneSurfNum];
337 0 : CarrollMRTDenominatorWin += zone_info.Fp[SendZoneSurfNum] * zone_info.Area[SendZoneSurfNum];
338 : }
339 : }
340 0 : if (CarrollMRTDenominatorWin > 0.0) {
341 0 : CarrollMRTInKTo4thWin = pow_4(CarrollMRTNumeratorWin / CarrollMRTDenominatorWin + Constant::Kelvin);
342 : }
343 0 : state.dataSurface->SurfWinIRfromParentZone(RecSurfNum) +=
344 0 : (zone_info.Fp[RecZoneSurfNum] * CarrollMRTInKTo4thWin) / SurfaceEmiss[RecZoneSurfNum];
345 : }
346 0 : netLWRadToRecSurf += zone_info.Fp[RecZoneSurfNum] * (CarrollMRTInKTo4th - SurfaceTempInKto4th[RecZoneSurfNum]);
347 : }
348 : } else {
349 13788155 : for (size_type RecZoneSurfNum = 0; RecZoneSurfNum < s_zone_Surfaces; ++RecZoneSurfNum) {
350 11795833 : int const RecSurfNum = zone_info.SurfacePtr[RecZoneSurfNum];
351 11795833 : int const ConstrNumRec = state.dataSurface->Surface(RecSurfNum).Construction;
352 11795833 : auto const &rec_construct = state.dataConstruction->Construct(ConstrNumRec);
353 11795833 : auto &netLWRadToRecSurf = NetLWRadToSurf(RecSurfNum);
354 :
355 : // Calculate net long-wave radiation for opaque surfaces and incident
356 : // long-wave radiation for windows.
357 11795833 : if (rec_construct.TypeIsWindow) { // Window
358 : // auto& rec_surface_window(SurfaceWindow(RecSurfNum));
359 549221 : Real64 scriptF_acc(0.0); // Local accumulator
360 549221 : Real64 netLWRadToRecSurf_cor(0.0); // Correction
361 549221 : Real64 IRfromParentZone_acc(0.0); // Local accumulator
362 7181532 : for (size_type SendZoneSurfNum = 0; SendZoneSurfNum < s_zone_Surfaces; ++SendZoneSurfNum) {
363 6632311 : size_type lSR = RecZoneSurfNum * s_zone_Surfaces + SendZoneSurfNum;
364 6632311 : Real64 const scriptF(zone_ScriptF[lSR]); // [ lSR ] == ( SendZoneSurfNum+1, RecZoneSurfNum+1 )
365 6632311 : Real64 const scriptF_temp_ink_4th(scriptF * SurfaceTempInKto4th[SendZoneSurfNum]);
366 : // Calculate interior LW incident on window rather than net LW for use in window layer heat balance calculation.
367 6632311 : IRfromParentZone_acc += scriptF_temp_ink_4th;
368 :
369 6632311 : if (RecZoneSurfNum != SendZoneSurfNum) {
370 6083090 : scriptF_acc += scriptF;
371 : } else {
372 549221 : netLWRadToRecSurf_cor = scriptF_temp_ink_4th;
373 : }
374 :
375 : // Per BG -- this should never happened. (CR6346,CR6550 caused this to be put in. Now removed. LKL 1/2013)
376 : // IF (SurfaceWindow(RecSurfNum)%IRfromParentZone < 0.0) THEN
377 : // CALL ShowRecurringWarningErrorAtEnd(state, 'CalcInteriorRadExchange: Window_IRFromParentZone negative,
378 : // Window="'// &
379 : // TRIM(Surface(RecSurfNum)%Name)//'"', &
380 : // SurfaceWindow(RecSurfNum)%IRErrCount)
381 : // CALL ShowRecurringContinueErrorAtEnd(state, '..occurs in Zone="'//TRIM(Surface(RecSurfNum)%ZoneName)// &
382 : // '", reset to 0.0 for remaining calculations.',SurfaceWindow(RecSurfNum)%IRErrCountC)
383 : // SurfaceWindow(RecSurfNum)%IRfromParentZone=0.0
384 : // ENDIF
385 : }
386 549221 : netLWRadToRecSurf += IRfromParentZone_acc - netLWRadToRecSurf_cor - (scriptF_acc * SurfaceTempInKto4th[RecZoneSurfNum]);
387 549221 : state.dataSurface->SurfWinIRfromParentZone(RecSurfNum) += IRfromParentZone_acc / SurfaceEmiss[RecZoneSurfNum];
388 : } else {
389 11246612 : Real64 netLWRadToRecSurf_acc(0.0); // Local accumulator
390 11246612 : zone_ScriptF[RecZoneSurfNum * s_zone_Surfaces + RecZoneSurfNum] = 0;
391 85870380 : for (size_type SendZoneSurfNum = 0; SendZoneSurfNum < s_zone_Surfaces; ++SendZoneSurfNum) {
392 74623768 : size_type lSR = RecZoneSurfNum * s_zone_Surfaces + SendZoneSurfNum;
393 74623768 : netLWRadToRecSurf_acc +=
394 149247536 : zone_ScriptF[lSR] * (SurfaceTempInKto4th[SendZoneSurfNum] -
395 74623768 : SurfaceTempInKto4th[RecZoneSurfNum]); // [ lSR ] == ( SendZoneSurfNum+1, RecZoneSurfNum+1 )
396 : }
397 11246612 : netLWRadToRecSurf += netLWRadToRecSurf_acc;
398 : }
399 : }
400 : }
401 : }
402 :
403 : // Automatic Surface Multipliers: Update values of surfaces not simulated
404 1362285 : if (state.dataSurface->UseRepresentativeSurfaceCalculations) {
405 0 : for (int SurfNum : state.dataSurface->AllHTSurfaceList) {
406 0 : int RepSurfNum = state.dataSurface->Surface(SurfNum).RepresentativeCalcSurfNum;
407 0 : if (SurfNum != RepSurfNum) {
408 0 : state.dataSurface->SurfWinIRfromParentZone(SurfNum) = state.dataSurface->SurfWinIRfromParentZone(RepSurfNum);
409 0 : NetLWRadToSurf(SurfNum) = NetLWRadToSurf(RepSurfNum);
410 : }
411 : }
412 : }
413 : }
414 :
415 3 : void UpdateMovableInsulationFlag(EnergyPlusData &state, bool &change, int const SurfNum)
416 : {
417 :
418 : // SUBROUTINE INFORMATION:
419 : // AUTHOR Rick Strand
420 : // DATE WRITTEN July 2016
421 :
422 : // PURPOSE OF THIS SUBROUTINE:
423 : // To determine if any changes in interior movable insulation have happened.
424 : // If there have been changes due to a schedule change AND a change in properties,
425 : // then the matrices which are used to calculate interior radiation must be recalculated.
426 :
427 3 : auto &s_surf = state.dataSurface;
428 :
429 3 : change = false;
430 :
431 3 : auto &movInsul = s_surf->intMovInsuls(SurfNum);
432 3 : if (movInsul.present != movInsul.presentPrevTS)
433 2 : change = (std::abs(state.dataConstruction->Construct(s_surf->Surface(SurfNum).Construction).InsideAbsorpThermal -
434 2 : state.dataMaterial->materials(movInsul.matNum)->AbsorpThermal) > 0.01);
435 3 : }
436 :
437 113 : void InitInteriorRadExchange(EnergyPlusData &state)
438 : {
439 :
440 : // SUBROUTINE INFORMATION:
441 : // AUTHOR Rick Strand
442 : // DATE WRITTEN September 2000
443 :
444 : // PURPOSE OF THIS SUBROUTINE:
445 : // Initializes the various parameters for Hottel's ScriptF method for
446 : // the grey interchange between surfaces in an enclosure.
447 :
448 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
449 : static constexpr std::string_view RoutineName("InitInteriorRadExchange: ");
450 : bool NoUserInputF; // Logical flag signifying no input F's for zone
451 113 : bool ErrorsFound(false);
452 : Real64 CheckValue1;
453 : Real64 CheckValue2;
454 : Real64 FinalCheckValue;
455 113 : Array2D<Real64> SaveApproximateViewFactors; // Save for View Factor reporting
456 : Real64 RowSum;
457 : Real64 FixedRowSum;
458 : int NumIterations;
459 113 : std::string Option1; // view factor report option
460 :
461 113 : auto &ViewFactorReport = state.dataHeatBalIntRadExchg->ViewFactorReport;
462 :
463 339 : General::ScanForReports(state, "ViewFactorInfo", ViewFactorReport, _, Option1);
464 :
465 113 : if (ViewFactorReport) { // Print heading
466 0 : print(state.files.eio, "{}\n", "! <Surface View Factor and Grey Interchange Information>");
467 0 : print(state.files.eio, "{}\n", "! <View Factor - Zone/Enclosure Information>,Zone/Enclosure Name,Number of Surfaces");
468 0 : print(state.files.eio,
469 : "{}\n",
470 : "! <View Factor - Surface Information>,Surface Name,Surface Class,Area {m2},Azimuth,Tilt,Thermal Emissivity,#Sides,Vertices");
471 0 : print(state.files.eio, "{}\n", "! <View Factor / Grey Interchange Type>,Surface Name(s)");
472 0 : print(state.files.eio, "{}\n", "! <View Factor>,Surface Name,Surface Class,Row Sum,View Factors for each Surface");
473 : }
474 :
475 113 : state.dataHeatBalIntRadExchg->MaxNumOfRadEnclosureSurfs = 0;
476 256 : for (int enclosureNum = 1; enclosureNum <= state.dataViewFactor->NumOfRadiantEnclosures; ++enclosureNum) {
477 143 : auto &thisEnclosure(state.dataViewFactor->EnclRadInfo(enclosureNum));
478 143 : if (enclosureNum == 1) {
479 100 : if (state.dataGlobal->DisplayAdvancedReportVariables) {
480 0 : print(state.files.eio,
481 : "{}\n",
482 : "! <Surface View Factor Check Values>,Zone/Enclosure Name,Original Check Value,Calculated Fixed Check Value,Final Check "
483 : "Value,Number of Iterations,Fixed RowSum Convergence,Used RowSum Convergence");
484 : }
485 : }
486 143 : int numEnclosureSurfaces = 0;
487 296 : for (int spaceNum : thisEnclosure.spaceNums) {
488 : // Note that Space.Surfaces only includes HT surfs, see SurfaceGeometry::CreateMissingSpaces
489 : // But it also includes air boundary surfaces which need to be excluded here
490 954 : for (int surfNum : state.dataHeatBal->space(spaceNum).surfaces) {
491 801 : if (state.dataSurface->Surface(surfNum).IsAirBoundarySurf) continue;
492 : // Automatic Surface Multipliers: Only include representative surfaces
493 799 : if (surfNum == state.dataSurface->Surface(surfNum).RepresentativeCalcSurfNum) {
494 799 : ++numEnclosureSurfaces;
495 : }
496 : }
497 : }
498 143 : thisEnclosure.NumOfSurfaces = numEnclosureSurfaces;
499 143 : state.dataHeatBalIntRadExchg->MaxNumOfRadEnclosureSurfs =
500 143 : max(state.dataHeatBalIntRadExchg->MaxNumOfRadEnclosureSurfs, numEnclosureSurfaces);
501 143 : if (numEnclosureSurfaces < 1) {
502 0 : ShowSevereError(state, format("{}No surfaces in enclosure={}.", RoutineName, thisEnclosure.Name));
503 0 : ErrorsFound = true;
504 : }
505 :
506 : // Allocate the parts of the derived type
507 143 : thisEnclosure.F.dimension(numEnclosureSurfaces, numEnclosureSurfaces, 0.0);
508 143 : thisEnclosure.ScriptF.dimension(numEnclosureSurfaces, numEnclosureSurfaces, 0.0);
509 143 : thisEnclosure.Area.dimension(numEnclosureSurfaces, 0.0);
510 143 : thisEnclosure.Emissivity.dimension(numEnclosureSurfaces, 0.0);
511 143 : thisEnclosure.Azimuth.dimension(numEnclosureSurfaces, 0.0);
512 143 : thisEnclosure.Tilt.dimension(numEnclosureSurfaces, 0.0);
513 143 : thisEnclosure.Fp.dimension(numEnclosureSurfaces, 1.0);
514 143 : thisEnclosure.FMRT.dimension(numEnclosureSurfaces, 0.0);
515 143 : thisEnclosure.SurfacePtr.dimension(numEnclosureSurfaces, 0);
516 :
517 : // Initialize the enclosure surface arrays
518 143 : int enclosureSurfNum = 0;
519 296 : for (int const spaceNum : thisEnclosure.spaceNums) {
520 153 : int priorZoneTotEnclSurfs = enclosureSurfNum;
521 954 : for (int surfNum : state.dataHeatBal->space(spaceNum).surfaces) {
522 801 : if (state.dataSurface->Surface(surfNum).IsAirBoundarySurf) continue;
523 : // Automatic Surface Multipliers: Only include representative surfaces
524 799 : if (surfNum == state.dataSurface->Surface(surfNum).RepresentativeCalcSurfNum) {
525 799 : ++enclosureSurfNum;
526 799 : thisEnclosure.SurfacePtr(enclosureSurfNum) = surfNum;
527 799 : thisEnclosure.Area(enclosureSurfNum) = state.dataSurface->Surface(surfNum).Area;
528 1598 : thisEnclosure.Emissivity(enclosureSurfNum) =
529 799 : state.dataConstruction->Construct(state.dataSurface->Surface(surfNum).Construction).InsideAbsorpThermal;
530 799 : thisEnclosure.Azimuth(enclosureSurfNum) = state.dataSurface->Surface(surfNum).Azimuth;
531 799 : thisEnclosure.Tilt(enclosureSurfNum) = state.dataSurface->Surface(surfNum).Tilt;
532 : }
533 : }
534 :
535 954 : for (int surfNum : state.dataHeatBal->space(spaceNum).surfaces) {
536 801 : if (state.dataSurface->Surface(surfNum).IsAirBoundarySurf) continue;
537 : // Map non-representative surfaces
538 799 : if (surfNum != state.dataSurface->Surface(surfNum).RepresentativeCalcSurfNum) {
539 : // Automatic Surface Multipliers: search for corresponding enclosure surface number
540 0 : for (int enclSNum = priorZoneTotEnclSurfs + 1; enclSNum <= enclosureSurfNum; ++enclSNum) {
541 0 : if (thisEnclosure.SurfacePtr(enclSNum) == state.dataSurface->Surface(surfNum).RepresentativeCalcSurfNum) {
542 : // enclosureSurfMap[allSurfNum] = enclSNum;
543 : // Increase the area for the combined enclosure surface
544 0 : thisEnclosure.Area(enclSNum) += state.dataSurface->Surface(surfNum).Area;
545 : }
546 : }
547 : }
548 : // Store SurfaceReportNums to maintain original reporting order
549 3020 : for (int enclSNum = priorZoneTotEnclSurfs + 1; enclSNum <= enclosureSurfNum; ++enclSNum) {
550 3020 : if (thisEnclosure.SurfacePtr(enclSNum) == state.dataSurface->AllSurfaceListReportOrder[surfNum - 1]) {
551 799 : thisEnclosure.SurfaceReportNums.push_back(enclSNum);
552 799 : break;
553 : }
554 : }
555 : }
556 : }
557 :
558 143 : if (thisEnclosure.NumOfSurfaces == 1) {
559 : // If there is only one surface in a zone, then there is no radiant exchange
560 5 : thisEnclosure.F = 0.0;
561 5 : thisEnclosure.ScriptF = 0.0;
562 5 : thisEnclosure.Fp = 0.0;
563 5 : thisEnclosure.FMRT = 0.0;
564 5 : if (state.dataGlobal->DisplayAdvancedReportVariables)
565 0 : print(state.files.eio, "Surface View Factor Check Values,{},0,0,0,-1,0,0\n", thisEnclosure.Name);
566 5 : continue; // Go to the next enclosure in the loop
567 : }
568 :
569 : // Process View Factors
570 138 : if (state.dataHeatBalIntRadExchg->CarrollMethod) {
571 :
572 : // User View Factors cannot be used with Carroll method.
573 0 : if (state.dataInputProcessing->inputProcessor->getNumObjectsFound(state, "ZoneProperty:UserViewFactors:BySurfaceName")) {
574 0 : ShowWarningError(state, "ZoneProperty:UserViewFactors:BySurfaceName objects have been defined, however View");
575 0 : ShowContinueError(state, " Factors are not used when Zone Radiant Exchange Algorithm is set to CarrollMRT.");
576 : }
577 0 : CalcFMRT(state, thisEnclosure.NumOfSurfaces, thisEnclosure.Area, thisEnclosure.FMRT);
578 0 : CalcFp(thisEnclosure.NumOfSurfaces, thisEnclosure.Emissivity, thisEnclosure.FMRT, thisEnclosure.Fp);
579 : } else {
580 : // Get user supplied view factors if available in idf.
581 :
582 138 : NoUserInputF = true;
583 :
584 138 : constexpr std::string_view cCurrentModuleObject = "ZoneProperty:UserViewFactors:BySurfaceName";
585 138 : int NumZonesWithUserFbyS = state.dataInputProcessing->inputProcessor->getNumObjectsFound(state, cCurrentModuleObject);
586 138 : if (NumZonesWithUserFbyS > 0) {
587 :
588 0 : GetInputViewFactorsbyName(state,
589 0 : thisEnclosure.Name,
590 : thisEnclosure.NumOfSurfaces,
591 : thisEnclosure.F,
592 0 : thisEnclosure.SurfacePtr,
593 : NoUserInputF,
594 : ErrorsFound); // Obtains user input view factors from input file
595 : }
596 :
597 138 : if (NoUserInputF) {
598 :
599 : // Calculate the view factors and make sure they satisfy reciprocity
600 138 : CalcApproximateViewFactors(state,
601 : thisEnclosure.NumOfSurfaces,
602 138 : thisEnclosure.Area,
603 138 : thisEnclosure.Azimuth,
604 138 : thisEnclosure.Tilt,
605 : thisEnclosure.F,
606 138 : thisEnclosure.SurfacePtr);
607 : }
608 :
609 138 : if (ViewFactorReport) { // Allocate and save user or approximate view factors for reporting.
610 0 : SaveApproximateViewFactors.allocate(thisEnclosure.NumOfSurfaces, thisEnclosure.NumOfSurfaces);
611 0 : SaveApproximateViewFactors = thisEnclosure.F;
612 : }
613 :
614 138 : bool anyIntMassInZone = DoesZoneHaveInternalMass(state, thisEnclosure.NumOfSurfaces, thisEnclosure.SurfacePtr);
615 138 : FixViewFactors(state,
616 : thisEnclosure.NumOfSurfaces,
617 138 : thisEnclosure.Area,
618 : thisEnclosure.F,
619 138 : thisEnclosure.Name,
620 138 : thisEnclosure.spaceNums,
621 : CheckValue1,
622 : CheckValue2,
623 : FinalCheckValue,
624 : NumIterations,
625 : FixedRowSum,
626 : anyIntMassInZone);
627 :
628 : // Calculate the script F factors
629 138 : CalcScriptF(state, thisEnclosure.NumOfSurfaces, thisEnclosure.Area, thisEnclosure.F, thisEnclosure.Emissivity, thisEnclosure.ScriptF);
630 138 : if (ViewFactorReport) { // Write to SurfInfo File
631 : // Zone Surface Information Output
632 0 : print(
633 0 : state.files.eio, "Surface View Factor - Zone/Enclosure Information,{},{}\n", thisEnclosure.Name, thisEnclosure.NumOfSurfaces);
634 :
635 0 : for (int SurfNum : thisEnclosure.SurfaceReportNums) {
636 0 : print(state.files.eio,
637 : "Surface View Factor - Surface Information,{},{},{:.4R},{:.4R},{:.4R},{:.4R},{}",
638 0 : state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Name,
639 0 : cSurfaceClass(state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Class),
640 : thisEnclosure.Area(SurfNum),
641 : thisEnclosure.Azimuth(SurfNum),
642 : thisEnclosure.Tilt(SurfNum),
643 : thisEnclosure.Emissivity(SurfNum),
644 0 : state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Sides);
645 0 : for (int Vindex = 1; Vindex <= state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Sides; ++Vindex) {
646 0 : auto const &Vertex = state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Vertex(Vindex);
647 0 : print(state.files.eio, ",{:.4R},{:.4R},{:.4R}", Vertex.x, Vertex.y, Vertex.z);
648 : }
649 0 : print(state.files.eio, "\n");
650 : }
651 :
652 0 : print(state.files.eio, "Approximate or User Input ViewFactors,To Surface,Surface Class,RowSum");
653 0 : for (int SurfNum : thisEnclosure.SurfaceReportNums) {
654 0 : print(state.files.eio, ",{}", state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Name);
655 : }
656 0 : print(state.files.eio, "\n");
657 :
658 0 : for (int Findex : thisEnclosure.SurfaceReportNums) {
659 0 : RowSum = sum(SaveApproximateViewFactors(_, Findex));
660 0 : print(state.files.eio,
661 : "{},{},{},{:.4R}",
662 : "View Factor",
663 0 : state.dataSurface->Surface(thisEnclosure.SurfacePtr(Findex)).Name,
664 0 : cSurfaceClass(state.dataSurface->Surface(thisEnclosure.SurfacePtr(Findex)).Class),
665 : RowSum);
666 0 : for (int SurfNum : thisEnclosure.SurfaceReportNums) {
667 0 : print(state.files.eio, ",{:.4R}", SaveApproximateViewFactors(SurfNum, Findex));
668 : }
669 0 : print(state.files.eio, "\n");
670 : }
671 :
672 0 : print(state.files.eio, "Final ViewFactors,To Surface,Surface Class,RowSum");
673 0 : for (int SurfNum : thisEnclosure.SurfaceReportNums) {
674 0 : print(state.files.eio, ",{}", state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Name);
675 : }
676 0 : print(state.files.eio, "\n");
677 :
678 0 : for (int Findex : thisEnclosure.SurfaceReportNums) {
679 0 : RowSum = sum(thisEnclosure.F(_, Findex));
680 0 : print(state.files.eio,
681 : "{},{},{},{:.4R}",
682 : "View Factor",
683 0 : state.dataSurface->Surface(thisEnclosure.SurfacePtr(Findex)).Name,
684 0 : cSurfaceClass(state.dataSurface->Surface(thisEnclosure.SurfacePtr(Findex)).Class),
685 : RowSum);
686 0 : for (int SurfNum : thisEnclosure.SurfaceReportNums) {
687 0 : print(state.files.eio, ",{:.4R}", thisEnclosure.F(SurfNum, Findex));
688 : }
689 0 : print(state.files.eio, "\n");
690 : }
691 :
692 0 : if (Option1 == "IDF") {
693 : // TODO Both "original" and "final" print the same output. This is likely a bug
694 : // (discovered while updating output to {fmt}
695 : // see:
696 : // https://github.com/NREL/EnergyPlusArchive/commit/1c08247853c297dce59f3f53cde47ccfa67720c0#diff-124964a7e9b73ce494c1952ab1acdeeb
697 0 : print(state.files.debug, "{}\n", "!======== original input factors ===========================");
698 0 : print(state.files.debug, "ZoneProperty:UserViewFactors:BySurfaceName,{},\n", thisEnclosure.Name);
699 0 : for (int SurfNum : thisEnclosure.SurfaceReportNums) {
700 0 : for (int Findex : thisEnclosure.SurfaceReportNums) {
701 0 : print(state.files.debug,
702 : " {},{},{:.6R}",
703 0 : state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Name,
704 0 : state.dataSurface->Surface(thisEnclosure.SurfacePtr(Findex)).Name,
705 : thisEnclosure.F(Findex, SurfNum));
706 0 : if (!(SurfNum == thisEnclosure.NumOfSurfaces && Findex == thisEnclosure.NumOfSurfaces)) {
707 0 : print(state.files.debug, ",\n");
708 : } else {
709 0 : print(state.files.debug, ";\n");
710 : }
711 : }
712 : }
713 0 : print(state.files.debug, "{}\n", "!============= end of data ======================");
714 :
715 0 : print(state.files.debug, "{}\n", "!============ final view factors =======================");
716 0 : print(state.files.debug, "ZoneProperty:UserViewFactors:BySurfaceName,{},\n", thisEnclosure.Name);
717 0 : for (int SurfNum : thisEnclosure.SurfaceReportNums) {
718 0 : for (int Findex : thisEnclosure.SurfaceReportNums) {
719 0 : print(state.files.debug,
720 : " {},{},{:.6R}",
721 0 : state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Name,
722 0 : state.dataSurface->Surface(thisEnclosure.SurfacePtr(Findex)).Name,
723 : thisEnclosure.F(Findex, SurfNum));
724 0 : if (!(SurfNum == thisEnclosure.SurfaceReportNums.back() && Findex == thisEnclosure.SurfaceReportNums.back())) {
725 0 : print(state.files.debug, ",\n");
726 : } else {
727 0 : print(state.files.debug, ";\n");
728 : }
729 : }
730 : }
731 0 : print(state.files.debug, "{}\n", "!============= end of data ======================");
732 : }
733 :
734 0 : print(state.files.eio, "Script F Factors,X Surface");
735 0 : for (int SurfNum : thisEnclosure.SurfaceReportNums) {
736 0 : print(state.files.eio, ",{}", state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Name);
737 : }
738 0 : print(state.files.eio, "\n");
739 0 : for (int Findex : thisEnclosure.SurfaceReportNums) {
740 0 : print(state.files.eio, "{},{}", "Script F Factor", state.dataSurface->Surface(thisEnclosure.SurfacePtr(Findex)).Name);
741 0 : for (int SurfNum : thisEnclosure.SurfaceReportNums) {
742 0 : print(state.files.eio, ",{:.4R}", thisEnclosure.ScriptF(Findex, SurfNum));
743 : }
744 0 : print(state.files.eio, "\n");
745 : }
746 :
747 : // Deallocate saved approximate/user view factors
748 0 : SaveApproximateViewFactors.deallocate();
749 : }
750 :
751 138 : RowSum = 0.0;
752 932 : for (int Findex : thisEnclosure.SurfaceReportNums) {
753 794 : RowSum += sum(thisEnclosure.F(_, Findex));
754 : }
755 138 : RowSum = std::abs(RowSum - thisEnclosure.NumOfSurfaces);
756 138 : FixedRowSum = std::abs(FixedRowSum - thisEnclosure.NumOfSurfaces);
757 138 : if (state.dataGlobal->DisplayAdvancedReportVariables) {
758 0 : print(state.files.eio,
759 : "Surface View Factor Check Values,{},{:.6R},{:.6R},{:.6R},{},{:.6R},{:.6R}\n",
760 0 : thisEnclosure.Name,
761 : CheckValue1,
762 : CheckValue2,
763 : FinalCheckValue,
764 : NumIterations,
765 : FixedRowSum,
766 : RowSum);
767 : }
768 : }
769 : }
770 :
771 113 : if (ErrorsFound) {
772 0 : ShowFatalError(state, format("{}Errors found during initialization of radiant exchange. Program terminated.", RoutineName));
773 : }
774 113 : }
775 :
776 135 : void InitSolarViewFactors(EnergyPlusData &state)
777 : {
778 :
779 : // Initializes view factors for diffuse solar distribution between surfaces in an enclosure.
780 :
781 135 : Array2D<Real64> SaveApproximateViewFactors; // Save for View Factor reporting
782 135 : std::string Option1; // view factor report option
783 : static constexpr std::string_view RoutineName("InitSolarViewFactors: ");
784 :
785 135 : bool ErrorsFound = false;
786 135 : bool ViewFactorReport = false;
787 405 : General::ScanForReports(state, "ViewFactorInfo", ViewFactorReport, _, Option1);
788 :
789 135 : if (ViewFactorReport) { // Print heading
790 0 : print(state.files.eio, "{}\n", "! <Solar View Factor Information>");
791 0 : print(state.files.eio, "{}\n", "! <Solar View Factor - Zone/Enclosure Information>,Zone/Enclosure Name,Number of Surfaces");
792 0 : print(state.files.eio,
793 : "{}\n",
794 : "! <Solar View Factor - Surface Information>,Surface Name,Surface Class,Area {m2},Azimuth,Tilt,Solar Absorbtance,#Sides,Vertices");
795 0 : print(state.files.eio, "{}\n", "! <Solar View Factor / Interchange Type>,Surface Name(s)");
796 0 : print(state.files.eio, "{}\n", "! <Solar View Factor>,Surface Name,Surface Class,Row Sum,View Factors for each Surface");
797 : }
798 :
799 135 : const std::string cCurrentModuleObject = "ZoneProperty:UserViewFactors:BySurfaceName";
800 135 : int NumZonesWithUserFbyS = state.dataInputProcessing->inputProcessor->getNumObjectsFound(state, cCurrentModuleObject);
801 135 : if (NumZonesWithUserFbyS > 0) AlignInputViewFactors(state, cCurrentModuleObject, ErrorsFound);
802 :
803 315 : for (int enclosureNum = 1; enclosureNum <= state.dataViewFactor->NumOfSolarEnclosures; ++enclosureNum) {
804 180 : auto &thisEnclosure = state.dataViewFactor->EnclSolInfo(enclosureNum);
805 180 : if (enclosureNum == 1) {
806 122 : if (state.dataGlobal->DisplayAdvancedReportVariables)
807 0 : print(state.files.eio,
808 : "{}\n",
809 : "! <Solar View Factor Check Values>,Zone/Enclosure Name,Original Check Value,Calculated Fixed Check "
810 : "Value,Final Check Value,Number of Iterations,Fixed RowSum Convergence,Used RowSum "
811 : "Convergence");
812 : }
813 180 : int numEnclosureSurfaces = 0;
814 371 : for (int spaceNum : thisEnclosure.spaceNums) {
815 : // Note that Space.Surfaces only includes HT surfs, see SurfaceGeometry::CreateMissingSpaces
816 : // But it also includes air boundary surfaces which need to be excluded here
817 1202 : for (int surfNum : state.dataHeatBal->space(spaceNum).surfaces) {
818 1011 : if (state.dataSurface->Surface(surfNum).IsAirBoundarySurf) continue;
819 1009 : ++numEnclosureSurfaces;
820 : }
821 : }
822 180 : thisEnclosure.NumOfSurfaces = numEnclosureSurfaces;
823 180 : if (numEnclosureSurfaces < 1) {
824 0 : ShowSevereError(state, format("{}No surfaces in enclosure={}.", RoutineName, thisEnclosure.Name));
825 0 : ErrorsFound = true;
826 : }
827 :
828 : // Allocate the parts of the derived type
829 180 : thisEnclosure.F.dimension(numEnclosureSurfaces, numEnclosureSurfaces, 0.0);
830 180 : thisEnclosure.Area.dimension(numEnclosureSurfaces, 0.0);
831 180 : thisEnclosure.SolAbsorptance.dimension(numEnclosureSurfaces, 0.0);
832 180 : thisEnclosure.Azimuth.dimension(numEnclosureSurfaces, 0.0);
833 180 : thisEnclosure.Tilt.dimension(numEnclosureSurfaces, 0.0);
834 180 : thisEnclosure.SurfacePtr.dimension(numEnclosureSurfaces, 0);
835 :
836 : // Initialize the surface pointer array
837 180 : int enclosureSurfNum = 0;
838 371 : for (int const spaceNum : thisEnclosure.spaceNums) {
839 191 : int priorZoneTotEnclSurfs = enclosureSurfNum;
840 1202 : for (int surfNum : state.dataHeatBal->space(spaceNum).surfaces) {
841 1011 : if (state.dataSurface->Surface(surfNum).IsAirBoundarySurf) continue;
842 1009 : ++enclosureSurfNum;
843 1009 : thisEnclosure.SurfacePtr(enclosureSurfNum) = surfNum;
844 : // Store pointers back to here
845 1009 : state.dataSurface->Surface(surfNum).SolarEnclSurfIndex = enclosureSurfNum;
846 1009 : state.dataSurface->Surface(surfNum).SolarEnclIndex = enclosureNum;
847 1009 : state.dataSurface->Surface(surfNum).RadEnclIndex = enclosureNum; // Radiant and Solar enclosures are parallel for now
848 1009 : if ((state.dataConstruction->Construct(state.dataSurface->Surface(surfNum).Construction).TypeIsWindow) &&
849 1009 : (state.dataSurface->Surface(surfNum).ExtBoundCond > 0) &&
850 0 : (state.dataSurface->Surface(surfNum).BaseSurf != surfNum)) { // Interzone window present
851 0 : thisEnclosure.HasInterZoneWindow = true;
852 : }
853 : }
854 : // Store SurfaceReportNums to maintain original reporting order
855 1202 : for (int surfNum : state.dataHeatBal->space(spaceNum).surfaces) {
856 3832 : for (int enclSNum = priorZoneTotEnclSurfs + 1; enclSNum <= enclosureSurfNum; ++enclSNum) {
857 3830 : if (thisEnclosure.SurfacePtr(enclSNum) == state.dataSurface->AllSurfaceListReportOrder[surfNum - 1]) {
858 1009 : thisEnclosure.SurfaceReportNums.push_back(enclSNum);
859 1009 : break;
860 : }
861 : }
862 : }
863 : }
864 : // Initialize the area and related arrays
865 1189 : for (int enclSurfNum = 1; enclSurfNum <= thisEnclosure.NumOfSurfaces; ++enclSurfNum) {
866 1009 : int const SurfNum = thisEnclosure.SurfacePtr(enclSurfNum);
867 1009 : thisEnclosure.Area(enclSurfNum) = state.dataSurface->Surface(SurfNum).Area;
868 2018 : thisEnclosure.SolAbsorptance(enclSurfNum) =
869 1009 : state.dataConstruction->Construct(state.dataSurface->Surface(SurfNum).Construction).InsideAbsorpSolar;
870 1009 : thisEnclosure.Azimuth(enclSurfNum) = state.dataSurface->Surface(SurfNum).Azimuth;
871 1009 : thisEnclosure.Tilt(enclSurfNum) = state.dataSurface->Surface(SurfNum).Tilt;
872 : }
873 :
874 180 : if (thisEnclosure.NumOfSurfaces == 1) {
875 : // If there is only one surface in a zone, then there is no solar distribution
876 10 : if (state.dataGlobal->DisplayAdvancedReportVariables)
877 0 : print(state.files.eio, "Solar View Factor Check Values,{},0,0,0,-1,0,0\n", thisEnclosure.Name);
878 10 : continue; // Go to the next enclosure in the loop
879 : }
880 :
881 : // Get user supplied view factors if available in idf.
882 :
883 170 : bool NoUserInputF = true;
884 :
885 170 : if (NumZonesWithUserFbyS > 0) {
886 :
887 4 : GetInputViewFactorsbyName(state,
888 4 : thisEnclosure.Name,
889 : thisEnclosure.NumOfSurfaces,
890 : thisEnclosure.F,
891 4 : thisEnclosure.SurfacePtr,
892 : NoUserInputF,
893 : ErrorsFound); // Obtains user input view factors from input file
894 : }
895 :
896 170 : if (NoUserInputF) {
897 :
898 : // Calculate the view factors and make sure they satisfy reciprocity
899 168 : CalcApproximateViewFactors(state,
900 : thisEnclosure.NumOfSurfaces,
901 168 : thisEnclosure.Area,
902 168 : thisEnclosure.Azimuth,
903 168 : thisEnclosure.Tilt,
904 : thisEnclosure.F,
905 168 : thisEnclosure.SurfacePtr);
906 : }
907 :
908 170 : if (ViewFactorReport) { // Allocate and save user or approximate view factors for reporting.
909 0 : SaveApproximateViewFactors.allocate(thisEnclosure.NumOfSurfaces, thisEnclosure.NumOfSurfaces);
910 0 : SaveApproximateViewFactors = thisEnclosure.F;
911 : }
912 :
913 170 : Real64 CheckValue1 = 0.0;
914 170 : Real64 CheckValue2 = 0.0;
915 170 : Real64 FinalCheckValue = 0.0;
916 170 : Real64 FixedRowSum = 0.0;
917 170 : int NumIterations = 0;
918 :
919 170 : bool anyIntMassInZone = DoesZoneHaveInternalMass(state, thisEnclosure.NumOfSurfaces, thisEnclosure.SurfacePtr);
920 170 : FixViewFactors(state,
921 : thisEnclosure.NumOfSurfaces,
922 170 : thisEnclosure.Area,
923 : thisEnclosure.F,
924 170 : thisEnclosure.Name,
925 170 : thisEnclosure.spaceNums,
926 : CheckValue1,
927 : CheckValue2,
928 : FinalCheckValue,
929 : NumIterations,
930 : FixedRowSum,
931 : anyIntMassInZone);
932 :
933 170 : if (ViewFactorReport) { // Write to SurfInfo File
934 : // Zone Surface Information Output
935 0 : print(state.files.eio, "Solar View Factor - Zone/Enclosure Information,{},{}\n", thisEnclosure.Name, thisEnclosure.NumOfSurfaces);
936 :
937 0 : for (int SurfNum : thisEnclosure.SurfaceReportNums) {
938 0 : print(state.files.eio,
939 : "Solar View Factor - Surface Information,{},{},{:.4R},{:.4R},{:.4R},{:.4R},{}",
940 0 : state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Name,
941 0 : cSurfaceClass(state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Class),
942 : thisEnclosure.Area(SurfNum),
943 : thisEnclosure.Azimuth(SurfNum),
944 : thisEnclosure.Tilt(SurfNum),
945 : thisEnclosure.SolAbsorptance(SurfNum),
946 0 : state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Sides);
947 :
948 0 : for (int Vindex = 1; Vindex <= state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Sides; ++Vindex) {
949 0 : auto const &Vertex = state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Vertex(Vindex);
950 0 : print(state.files.eio, ",{:.4R},{:.4R},{:.4R}", Vertex.x, Vertex.y, Vertex.z);
951 : }
952 0 : print(state.files.eio, "\n");
953 : }
954 :
955 0 : print(state.files.eio, "Approximate or User Input Solar ViewFactors,To Surface,Surface Class,RowSum");
956 0 : for (int SurfNum : thisEnclosure.SurfaceReportNums) {
957 0 : print(state.files.eio, ",{}", state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Name);
958 : }
959 0 : print(state.files.eio, "\n");
960 :
961 0 : for (int Findex : thisEnclosure.SurfaceReportNums) {
962 0 : Real64 RowSum = sum(SaveApproximateViewFactors(_, Findex));
963 0 : print(state.files.eio,
964 : "Solar View Factor,{},{},{:.4R}",
965 0 : state.dataSurface->Surface(thisEnclosure.SurfacePtr(Findex)).Name,
966 0 : cSurfaceClass(state.dataSurface->Surface(thisEnclosure.SurfacePtr(Findex)).Class),
967 : RowSum);
968 0 : for (int SurfNum : thisEnclosure.SurfaceReportNums) {
969 0 : print(state.files.eio, ",{:.4R}", SaveApproximateViewFactors(SurfNum, Findex));
970 : }
971 0 : print(state.files.eio, "\n");
972 : }
973 :
974 0 : print(state.files.eio, "Final Solar ViewFactors,To Surface,Surface Class,RowSum");
975 0 : for (int SurfNum : thisEnclosure.SurfaceReportNums) {
976 0 : print(state.files.eio, ",{}", state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Name);
977 : }
978 0 : print(state.files.eio, "\n");
979 :
980 0 : for (int Findex : thisEnclosure.SurfaceReportNums) {
981 0 : Real64 RowSum = sum(thisEnclosure.F(_, Findex));
982 0 : print(state.files.eio,
983 : "{},{},{},{:.4R}",
984 : "Solar View Factor",
985 0 : state.dataSurface->Surface(thisEnclosure.SurfacePtr(Findex)).Name,
986 0 : cSurfaceClass(state.dataSurface->Surface(thisEnclosure.SurfacePtr(Findex)).Class),
987 : RowSum);
988 0 : for (int SurfNum : thisEnclosure.SurfaceReportNums) {
989 0 : print(state.files.eio, ",{:.4R}", thisEnclosure.F(SurfNum, Findex));
990 : }
991 0 : print(state.files.eio, "\n");
992 : }
993 :
994 0 : if (Option1 == "IDF") {
995 : // TODO Both "original" and "final" print the same output. This is likely a bug
996 : // see:
997 : // https://github.com/NREL/EnergyPlusArchive/commit/1c08247853c297dce59f3f53cde47ccfa67720c0#diff-124964a7e9b73ce494c1952ab1acdeeb
998 0 : print(state.files.debug, "{}\n", "!======== original input factors ===========================");
999 0 : print(state.files.debug, "ZoneProperty:UserViewFactors:BySurfaceName,{},\n", thisEnclosure.Name);
1000 0 : for (int SurfNum : thisEnclosure.SurfaceReportNums) {
1001 0 : for (int Findex : thisEnclosure.SurfaceReportNums) {
1002 0 : print(state.files.debug,
1003 : " {},{},{:.6R}",
1004 0 : state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Name,
1005 0 : state.dataSurface->Surface(thisEnclosure.SurfacePtr(Findex)).Name,
1006 : thisEnclosure.F(Findex, SurfNum));
1007 0 : if (!(SurfNum == thisEnclosure.NumOfSurfaces && Findex == thisEnclosure.NumOfSurfaces)) {
1008 0 : print(state.files.debug, ",\n");
1009 : } else {
1010 0 : print(state.files.debug, ";\n");
1011 : }
1012 : }
1013 : }
1014 0 : print(state.files.debug, "{}\n", "!============= end of data ======================");
1015 :
1016 0 : print(state.files.debug, "{}\n", "!============ final view factors =======================");
1017 0 : print(state.files.debug, "ZoneProperty:UserViewFactors:BySurfaceName,{},\n", thisEnclosure.Name);
1018 0 : for (int SurfNum : thisEnclosure.SurfaceReportNums) {
1019 0 : for (int Findex : thisEnclosure.SurfaceReportNums) {
1020 0 : print(state.files.debug,
1021 : " {},{},{:.6R}",
1022 0 : state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Name,
1023 0 : state.dataSurface->Surface(thisEnclosure.SurfacePtr(Findex)).Name,
1024 : thisEnclosure.F(Findex, SurfNum));
1025 0 : if (!(SurfNum == thisEnclosure.NumOfSurfaces && Findex == thisEnclosure.NumOfSurfaces)) {
1026 0 : print(state.files.debug, ",\n");
1027 : } else {
1028 0 : print(state.files.debug, ";\n");
1029 : }
1030 : }
1031 : }
1032 0 : print(state.files.debug, "{}\n", "!============= end of data ======================");
1033 : }
1034 0 : SaveApproximateViewFactors.deallocate();
1035 : }
1036 :
1037 170 : Real64 RowSum = 0.0;
1038 1169 : for (int Findex : thisEnclosure.SurfaceReportNums) {
1039 999 : RowSum += sum(thisEnclosure.F(_, Findex));
1040 : }
1041 170 : RowSum = std::abs(RowSum - thisEnclosure.NumOfSurfaces);
1042 170 : FixedRowSum = std::abs(FixedRowSum - thisEnclosure.NumOfSurfaces);
1043 170 : if (state.dataGlobal->DisplayAdvancedReportVariables) {
1044 0 : print(state.files.eio,
1045 : "Solar View Factor Check Values,{},{:.6R},{:.6R},{:.6R},{},{:.6R},{:.6R}\n",
1046 0 : thisEnclosure.Name,
1047 : CheckValue1,
1048 : CheckValue2,
1049 : FinalCheckValue,
1050 : NumIterations,
1051 : FixedRowSum,
1052 : RowSum);
1053 : }
1054 : }
1055 :
1056 135 : if (ErrorsFound) {
1057 0 : ShowFatalError(state, format("{}Errors found during initialization of diffuse solar distribution. Program terminated.", RoutineName));
1058 : }
1059 135 : }
1060 :
1061 0 : void GetInputViewFactors(EnergyPlusData &state,
1062 : std::string const &ZoneName, // Needed to check for user input view factors.
1063 : int const N, // NUMBER OF SURFACES
1064 : Array2A<Real64> F, // USER INPUT DIRECT VIEW FACTOR MATRIX (N X N)
1065 : [[maybe_unused]] const Array1D_int &SPtr, // pointer to actual surface number
1066 : bool &NoUserInputF, // Flag signifying no input F's for this
1067 : bool &ErrorsFound // True when errors are found in number of fields vs max args
1068 : )
1069 : {
1070 :
1071 : // SUBROUTINE INFORMATION:
1072 : // AUTHOR Curt Pedersen
1073 : // DATE WRITTEN September 2005
1074 : // MODIFIED Linda Lawrie;September 2010
1075 :
1076 : // PURPOSE OF THIS SUBROUTINE:
1077 : // This routine gets the user view factor info.
1078 :
1079 : // Argument array dimensioning
1080 0 : F.dim(N, N);
1081 : // EP_SIZE_CHECK(SPtr, N);
1082 :
1083 0 : NoUserInputF = true;
1084 0 : int UserFZoneIndex = state.dataInputProcessing->inputProcessor->getObjectItemNum(state, "ZoneProperty:UserViewFactors", ZoneName);
1085 0 : if (UserFZoneIndex > 0) {
1086 : int NumAlphas;
1087 : int NumNums;
1088 : int IOStat;
1089 0 : NoUserInputF = false;
1090 :
1091 0 : state.dataInputProcessing->inputProcessor->getObjectItem(state,
1092 : "ZoneProperty:UserViewFactors",
1093 : UserFZoneIndex,
1094 0 : state.dataIPShortCut->cAlphaArgs,
1095 : NumAlphas,
1096 0 : state.dataIPShortCut->rNumericArgs,
1097 : NumNums,
1098 : IOStat,
1099 0 : state.dataIPShortCut->lNumericFieldBlanks,
1100 0 : state.dataIPShortCut->lAlphaFieldBlanks,
1101 0 : state.dataIPShortCut->cAlphaFieldNames,
1102 0 : state.dataIPShortCut->cNumericFieldNames);
1103 :
1104 0 : if (NumNums < 3 * pow_2(N)) {
1105 0 : std::string_view cCurrentModuleObject = "ZoneProperty:UserViewFactors";
1106 0 : ShowSevereError(state, format("GetInputViewFactors: {}=\"{}\", not enough values.", cCurrentModuleObject, ZoneName));
1107 0 : ShowContinueError(state, format("...Number of input values [{}] is less than the required number=[{}].", NumNums, 3 * pow_2(N)));
1108 0 : ErrorsFound = true;
1109 0 : NumNums = 0;
1110 : }
1111 0 : F = 0.0;
1112 0 : for (int index = 1; index <= NumNums; index += 3) {
1113 0 : int inx1 = state.dataIPShortCut->rNumericArgs(index);
1114 0 : int inx2 = state.dataIPShortCut->rNumericArgs(index + 1);
1115 0 : F(inx2, inx1) = state.dataIPShortCut->rNumericArgs(index + 2);
1116 : }
1117 : }
1118 0 : }
1119 :
1120 6 : void AlignInputViewFactors(EnergyPlusData &state,
1121 : std::string const &cCurrentModuleObject, // Object type
1122 : bool &ErrorsFound // True when errors are found
1123 : )
1124 : {
1125 6 : auto const instances = state.dataInputProcessing->inputProcessor->epJSON.find(cCurrentModuleObject);
1126 6 : auto &instancesValue = instances.value();
1127 20 : for (auto instance = instancesValue.begin(); instance != instancesValue.end(); ++instance) {
1128 14 : auto const &fields = instance.value();
1129 14 : std::string const thisSpaceOrSpaceListName = Util::makeUPPER(fields.at("zone_or_zonelist_or_space_or_spacelist_name").get<std::string>());
1130 : // do not mark object as used here - let GetInputViewFactorsbyName do that
1131 :
1132 : // Look for matching radiant enclosure name (solar enclosures have the same names)
1133 14 : bool enclMatchFound = false;
1134 48 : for (int enclosureNum = 1; enclosureNum <= state.dataViewFactor->NumOfRadiantEnclosures; ++enclosureNum) {
1135 40 : auto &thisEnclosure = state.dataViewFactor->EnclRadInfo(enclosureNum);
1136 40 : if (Util::SameString(thisSpaceOrSpaceListName, thisEnclosure.Name)) {
1137 : // View factor space name matches enclosure name
1138 6 : enclMatchFound = true;
1139 6 : break;
1140 : }
1141 : }
1142 14 : if (enclMatchFound) continue; // We're done with this instance
1143 : // Find matching SpaceList name
1144 8 : int spaceListNum = Util::FindItemInList(thisSpaceOrSpaceListName, state.dataHeatBal->spaceList);
1145 8 : int zoneListNum = 0;
1146 8 : int inputZoneNum = 0;
1147 8 : size_t listSize = 0;
1148 8 : if (spaceListNum > 0) {
1149 4 : listSize = int(state.dataHeatBal->spaceList(spaceListNum).spaces.size());
1150 : } else {
1151 : // Check Zone and ZoneList for backwards compatibility
1152 4 : zoneListNum = Util::FindItemInList(thisSpaceOrSpaceListName, state.dataHeatBal->ZoneList);
1153 4 : if (zoneListNum > 0) {
1154 0 : for (int const zoneNum : state.dataHeatBal->ZoneList(zoneListNum).Zone) {
1155 0 : listSize += state.dataHeatBal->Zone(zoneNum).spaceIndexes.size();
1156 : }
1157 : } else {
1158 4 : inputZoneNum = Util::FindItemInList(thisSpaceOrSpaceListName, state.dataHeatBal->Zone);
1159 4 : if (inputZoneNum > 0) {
1160 0 : listSize = state.dataHeatBal->Zone(inputZoneNum).spaceIndexes.size();
1161 : }
1162 : }
1163 : }
1164 :
1165 8 : if (listSize > 0) {
1166 : // Look for radiant enclosure with same list of spaces (solar enclosures are the same, so this is sufficient)
1167 10 : for (int enclosureNum = 1; enclosureNum <= state.dataViewFactor->NumOfRadiantEnclosures; ++enclosureNum) {
1168 8 : auto &thisEnclosure = state.dataViewFactor->EnclRadInfo(enclosureNum);
1169 8 : bool anySpaceNotFound = false;
1170 : // If the number of enclosure spaces is not the same as the number of spacelist space, go to the next enclosure
1171 8 : if (listSize != thisEnclosure.spaceNums.size()) continue;
1172 6 : if (spaceListNum > 0) {
1173 10 : for (int sListSpaceNum : state.dataHeatBal->spaceList(spaceListNum).spaces) {
1174 : // Search for matching spaces
1175 8 : bool thisSpaceFound = false;
1176 18 : for (int enclSpaceNum : thisEnclosure.spaceNums) {
1177 14 : if (enclSpaceNum == sListSpaceNum) {
1178 4 : thisSpaceFound = true;
1179 4 : break;
1180 : }
1181 : }
1182 8 : if (!thisSpaceFound) {
1183 4 : anySpaceNotFound = true;
1184 4 : break;
1185 : }
1186 : }
1187 0 : } else if (zoneListNum > 0) {
1188 0 : for (int zListZoneNum : state.dataHeatBal->ZoneList(zoneListNum).Zone) {
1189 0 : for (int spaceNum : state.dataHeatBal->Zone(zListZoneNum).spaceIndexes) {
1190 : // Search for matching spaces
1191 0 : bool thisSpaceFound = false;
1192 0 : for (int enclSpaceNum : thisEnclosure.spaceNums) {
1193 0 : if (enclSpaceNum == spaceNum) {
1194 0 : thisSpaceFound = true;
1195 0 : break;
1196 : }
1197 : }
1198 0 : if (!thisSpaceFound) {
1199 0 : anySpaceNotFound = true;
1200 0 : break;
1201 : }
1202 : }
1203 0 : if (anySpaceNotFound) {
1204 0 : break;
1205 : }
1206 : }
1207 0 : } else if (inputZoneNum > 0) {
1208 0 : for (int spaceNum : state.dataHeatBal->Zone(inputZoneNum).spaceIndexes) {
1209 : // Search for matching spaces
1210 0 : bool thisSpaceFound = false;
1211 0 : for (int enclSpaceNum : thisEnclosure.spaceNums) {
1212 0 : if (enclSpaceNum == spaceNum) {
1213 0 : thisSpaceFound = true;
1214 0 : break;
1215 : }
1216 : }
1217 0 : if (!thisSpaceFound) {
1218 0 : anySpaceNotFound = true;
1219 0 : break;
1220 : }
1221 : }
1222 : }
1223 6 : if (anySpaceNotFound) {
1224 4 : continue; // On to the next enclosure
1225 : } else {
1226 2 : enclMatchFound = true;
1227 : // If matching SpaceList or ZoneList or Zone found, set the radiant and solar enclosure names to match
1228 2 : thisEnclosure.Name = thisSpaceOrSpaceListName;
1229 2 : state.dataViewFactor->EnclSolInfo(enclosureNum).Name = thisSpaceOrSpaceListName;
1230 2 : break; // We're done with radiant enclosures
1231 : }
1232 : }
1233 : }
1234 8 : if (!enclMatchFound) {
1235 6 : if (spaceListNum > 0) {
1236 4 : ShowSevereError(state,
1237 4 : format("AlignInputViewFactors: {}=\"{}\" found a matching SpaceList, but did not find a matching radiant or "
1238 : "solar enclosure with the same spaces.",
1239 : cCurrentModuleObject,
1240 : thisSpaceOrSpaceListName));
1241 2 : ErrorsFound = true;
1242 :
1243 4 : } else if (zoneListNum > 0) {
1244 0 : ShowSevereError(state,
1245 0 : format("AlignInputViewFactors: {}=\"{}\" found a matching ZoneList, but did not find a matching radiant or solar "
1246 : "enclosure with the same spaces.",
1247 : cCurrentModuleObject,
1248 : thisSpaceOrSpaceListName));
1249 0 : ErrorsFound = true;
1250 :
1251 : } else {
1252 8 : ShowSevereError(state,
1253 8 : format("AlignInputViewFactors: {}=\"{}\" did not find a matching radiant or solar enclosure name.",
1254 : cCurrentModuleObject,
1255 : thisSpaceOrSpaceListName));
1256 4 : ErrorsFound = true;
1257 : }
1258 : }
1259 14 : }
1260 6 : }
1261 :
1262 4 : void GetInputViewFactorsbyName(EnergyPlusData &state,
1263 : std::string const &EnclosureName, // Needed to check for user input view factors.
1264 : int const N, // NUMBER OF SURFACES
1265 : Array2A<Real64> F, // USER INPUT DIRECT VIEW FACTOR MATRIX (N X N)
1266 : const Array1D_int &SPtr, // pointer to actual surface number
1267 : bool &NoUserInputF, // Flag signifying no input F's for this
1268 : bool &ErrorsFound // True when errors are found in number of fields vs max args
1269 : )
1270 : {
1271 :
1272 : // SUBROUTINE INFORMATION:
1273 : // AUTHOR Curt Pedersen
1274 : // DATE WRITTEN September 2005
1275 : // MODIFIED Linda Lawrie;September 2010
1276 :
1277 : // PURPOSE OF THIS SUBROUTINE:
1278 : // This routine gets the user view factor info for an enclosure which could be a zone or a group of zones
1279 :
1280 : // Argument array dimensioning
1281 4 : F.dim(N, N);
1282 4 : EP_SIZE_CHECK(SPtr, N);
1283 :
1284 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
1285 : int UserFZoneIndex;
1286 4 : Array1D_string enclosureSurfaceNames;
1287 4 : NoUserInputF = true;
1288 8 : UserFZoneIndex = state.dataInputProcessing->inputProcessor->getObjectItemNum(
1289 : state, "ZoneProperty:UserViewFactors:BySurfaceName", "zone_or_zonelist_or_space_or_spacelist_name", EnclosureName);
1290 :
1291 4 : if (UserFZoneIndex > 0) {
1292 : int NumAlphas;
1293 : int NumNums;
1294 : int IOStat;
1295 2 : auto &cCurrentModuleObject = state.dataIPShortCut->cCurrentModuleObject;
1296 2 : enclosureSurfaceNames.allocate(N);
1297 16 : for (int index = 1; index <= N; ++index) {
1298 14 : enclosureSurfaceNames(index) = state.dataSurface->Surface(SPtr(index)).Name;
1299 : }
1300 2 : NoUserInputF = false;
1301 :
1302 4 : state.dataInputProcessing->inputProcessor->getObjectItem(state,
1303 : "ZoneProperty:UserViewFactors:BySurfaceName",
1304 : UserFZoneIndex,
1305 2 : state.dataIPShortCut->cAlphaArgs,
1306 : NumAlphas,
1307 2 : state.dataIPShortCut->rNumericArgs,
1308 : NumNums,
1309 : IOStat,
1310 2 : state.dataIPShortCut->lNumericFieldBlanks,
1311 2 : state.dataIPShortCut->lAlphaFieldBlanks,
1312 2 : state.dataIPShortCut->cAlphaFieldNames,
1313 2 : state.dataIPShortCut->cNumericFieldNames);
1314 :
1315 2 : F = 0.0;
1316 2 : int numinx1 = 0;
1317 2 : if (NumNums < pow_2(N)) {
1318 0 : ShowWarningError(state, format("GetInputViewFactors: {}=\"{}\", not enough values.", cCurrentModuleObject, EnclosureName));
1319 0 : ShowContinueError(state,
1320 0 : format("...Number of input values [{}] is less than the required number=[{}] Missing surface pairs will have a "
1321 : "zero view factor.",
1322 : NumNums,
1323 0 : pow_2(N)));
1324 : }
1325 :
1326 100 : for (int index = 2; index <= NumAlphas; index += 2) {
1327 98 : int inx1 = Util::FindItemInList(state.dataIPShortCut->cAlphaArgs(index), enclosureSurfaceNames, N);
1328 98 : int inx2 = Util::FindItemInList(state.dataIPShortCut->cAlphaArgs(index + 1), enclosureSurfaceNames, N);
1329 98 : if (inx1 == 0) {
1330 0 : ShowSevereError(state, format("GetInputViewFactors: {}=\"{}\", invalid surface name.", cCurrentModuleObject, EnclosureName));
1331 0 : ShowContinueError(state,
1332 0 : format("...Surface name=\"{}\", not in this zone or enclosure.", state.dataIPShortCut->cAlphaArgs(index)));
1333 0 : ErrorsFound = true;
1334 : }
1335 98 : if (inx2 == 0) {
1336 0 : ShowSevereError(state, format("GetInputViewFactors: {}=\"{}\", invalid surface name.", cCurrentModuleObject, EnclosureName));
1337 0 : ShowContinueError(state,
1338 0 : format("...Surface name=\"{}\", not in this zone or enclosure.", state.dataIPShortCut->cAlphaArgs(index + 2)));
1339 0 : ErrorsFound = true;
1340 : }
1341 98 : ++numinx1;
1342 98 : if (inx1 > 0 && inx2 > 0) F(inx2, inx1) = state.dataIPShortCut->rNumericArgs(numinx1);
1343 : }
1344 :
1345 2 : enclosureSurfaceNames.deallocate();
1346 : }
1347 4 : }
1348 :
1349 307 : void CalcApproximateViewFactors(EnergyPlusData &state,
1350 : int const N, // NUMBER OF SURFACES
1351 : const Array1D<Real64> &A, // AREA VECTOR- ASSUMED,BE N ELEMENTS LONG
1352 : const Array1D<Real64> &Azimuth, // Facing angle of the surface (in degrees)
1353 : const Array1D<Real64> &Tilt, // Tilt angle of the surface (in degrees)
1354 : Array2A<Real64> F, // APPROXIMATE DIRECT VIEW FACTOR MATRIX (N X N)
1355 : const Array1D_int &SPtr // pointer to REAL(r64) surface number (for error message)
1356 : )
1357 : {
1358 :
1359 : // SUBROUTINE INFORMATION:
1360 : // AUTHOR Curt Pedersen
1361 : // DATE WRITTEN July 2000
1362 : // MODIFIED March 2001 (RKS) to disallow surfaces facing the same direction to interact radiatively
1363 : // May 2002 (COP) to include INTMASS, FLOOR, ROOF and CEILING.
1364 : // RE-ENGINEERED September 2000 (RKS for EnergyPlus)
1365 :
1366 : // PURPOSE OF THIS SUBROUTINE:
1367 : // This subroutine approximates view factors using an area weighting.
1368 : // This is improved by one degree by not allowing surfaces facing the same
1369 : // direction to "see" each other.
1370 :
1371 : // METHODOLOGY EMPLOYED:
1372 : // Each surface sees some area of other surfaces within the zone. The view
1373 : // factors from the surface to the other seen surfaces are defined by their
1374 : // area over the summed area of seen surfaces. Surfaces facing the same angle
1375 : // are assumed to not be able to see each other.
1376 : // Modified May 2002 to cover poorly defined surface orientation. Now all thermal masses, roofs and
1377 : // ceilings are "seen" by other surfaces. Floors are seen by all other surfaces, but
1378 : // not by other floors.
1379 :
1380 : // Argument array dimensioning
1381 307 : EP_SIZE_CHECK(A, N);
1382 307 : EP_SIZE_CHECK(Azimuth, N);
1383 307 : EP_SIZE_CHECK(Tilt, N);
1384 307 : F.dim(N, N);
1385 307 : EP_SIZE_CHECK(SPtr, N);
1386 :
1387 : // SUBROUTINE PARAMETER DEFINITIONS:
1388 307 : Real64 constexpr SameAngleLimit(10.0); // If the difference in the azimuth angles are above this value (degrees),
1389 : // then the surfaces are assumed to be facing different directions.
1390 :
1391 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
1392 : int i; // DO loop counters for surfaces in the zone
1393 : int j;
1394 307 : Array1D<Real64> ZoneArea; // Sum of the area of all zone surfaces seen
1395 :
1396 : // Calculate the sum of the areas seen by all zone surfaces
1397 307 : ZoneArea.dimension(N, 0.0);
1398 2095 : for (i = 1; i <= N; ++i) {
1399 13688 : for (j = 1; j <= N; ++j) {
1400 : // No surface sees itself and no floor sees another floor
1401 13630 : if (i == j || (state.dataSurface->Surface(SPtr(j)).Class == DataSurfaces::SurfaceClass::Floor &&
1402 1730 : state.dataSurface->Surface(SPtr(i)).Class == DataSurfaces::SurfaceClass::Floor))
1403 1938 : continue;
1404 : // All surface types see internal mass
1405 : // All surface types see floors
1406 : // Floors always see ceilings/roofs
1407 : // All other surfaces whose tilt or facing angle differences are greater than 10 degrees see each other
1408 9962 : if ((state.dataSurface->Surface(SPtr(j)).Class == DataSurfaces::SurfaceClass::IntMass) ||
1409 9890 : (state.dataSurface->Surface(SPtr(i)).Class == DataSurfaces::SurfaceClass::IntMass) ||
1410 9818 : (state.dataSurface->Surface(SPtr(j)).Class == DataSurfaces::SurfaceClass::Floor) ||
1411 14926 : (state.dataSurface->Surface(SPtr(i)).Class == DataSurfaces::SurfaceClass::Floor) ||
1412 26536 : (std::abs(Azimuth(i) - Azimuth(j)) > SameAngleLimit && std::abs(Azimuth(i) - Azimuth(j)) < 360.0 - SameAngleLimit) ||
1413 1018 : (std::abs(Tilt(i) - Tilt(j)) > SameAngleLimit)) {
1414 9598 : ZoneArea(i) += A(j);
1415 : }
1416 : }
1417 1788 : if (ZoneArea(i) <= 0.0) {
1418 76 : ShowWarningError(state, "CalcApproximateViewFactors: Zero area for all other zone surfaces.");
1419 76 : ShowContinueError(state,
1420 76 : format("Happens for Surface=\"{}\" in Zone={}",
1421 38 : state.dataSurface->Surface(SPtr(i)).Name,
1422 38 : state.dataHeatBal->Zone(state.dataSurface->Surface(SPtr(i)).Zone).Name));
1423 : }
1424 : }
1425 :
1426 : // Set up the approximate view factors. First these are initialized to all zero.
1427 : // This will clear out any junk leftover from whenever. Then, for each zone
1428 : // surface, set the view factor from that surface to other surfaces as the
1429 : // area of the other surface divided by the sum of the area of all zone surfaces
1430 : // that the original surface can actually see (calculated above). This will
1431 : // allow that the sum of all view factors from the original surface to all other
1432 : // surfaces will equal unity. F(I,J)=0 if I=J or if the surfaces face the same
1433 : // direction.
1434 307 : F = 0.0;
1435 2095 : for (i = 1; i <= N; ++i) {
1436 13688 : for (j = 1; j <= N; ++j) {
1437 : // No surface sees itself and no floor sees another floor
1438 13630 : if (i == j || (state.dataSurface->Surface(SPtr(j)).Class == DataSurfaces::SurfaceClass::Floor &&
1439 1730 : state.dataSurface->Surface(SPtr(i)).Class == DataSurfaces::SurfaceClass::Floor))
1440 1938 : continue;
1441 : // All surface types see internal mass
1442 : // All surface types see floors
1443 : // Floors always see ceilings/roofs
1444 : // All other surfaces whose tilt or facing angle differences are greater than 10 degrees see each other
1445 9962 : if ((state.dataSurface->Surface(SPtr(j)).Class == DataSurfaces::SurfaceClass::IntMass) ||
1446 9890 : (state.dataSurface->Surface(SPtr(i)).Class == DataSurfaces::SurfaceClass::IntMass) ||
1447 9818 : (state.dataSurface->Surface(SPtr(j)).Class == DataSurfaces::SurfaceClass::Floor) ||
1448 14926 : (state.dataSurface->Surface(SPtr(i)).Class == DataSurfaces::SurfaceClass::Floor) ||
1449 26536 : (std::abs(Azimuth(i) - Azimuth(j)) > SameAngleLimit && std::abs(Azimuth(i) - Azimuth(j)) < 360.0 - SameAngleLimit) ||
1450 1018 : (std::abs(Tilt(i) - Tilt(j)) > SameAngleLimit)) {
1451 : // avoid a divide by zero if there are no other surfaces in the zone that can be seen
1452 9598 : if (ZoneArea(i) > 0.0) F(j, i) = A(j) / (ZoneArea(i));
1453 : }
1454 : }
1455 : }
1456 :
1457 307 : ZoneArea.deallocate();
1458 307 : }
1459 :
1460 313 : void FixViewFactors(EnergyPlusData &state,
1461 : int const N, // NUMBER OF SURFACES
1462 : const Array1D<Real64> &A, // AREA VECTOR- ASSUMED,BE N ELEMENTS LONG
1463 : Array2A<Real64> F, // APPROXIMATE DIRECT VIEW FACTOR MATRIX (N X N)
1464 : std::string &enclName, // Name of Enclosure being fixed
1465 : std::vector<int> const &spaceNums, // Zones which are part of this enclosure
1466 : Real64 &OriginalCheckValue, // check of SUM(F) - N
1467 : Real64 &FixedCheckValue, // check after fixed of SUM(F) - N
1468 : Real64 &FinalCheckValue, // the one to go with
1469 : int &NumIterations, // number of iterations to fixed
1470 : Real64 &RowSum, // RowSum of Fixed
1471 : bool const anyIntMassInZone // are there any internal mass surfaces in the zone
1472 : )
1473 : {
1474 :
1475 : // SUBROUTINE INFORMATION:
1476 : // AUTHOR Curt Pedersen
1477 : // DATE WRITTEN July 2000
1478 : // MODIFIED September 2000 (RKS for EnergyPlus)
1479 : // April 2005,COP added capability to handle a
1480 : // surface larger than sum of all others (nonenclosure)
1481 : // by using a Fii view factor for that surface. Process is
1482 : // now much more robust and stable.
1483 :
1484 : // PURPOSE OF THIS SUBROUTINE:
1485 : // This subroutine fixes approximate view factors and enforces reciprocity
1486 : // and completeness.
1487 :
1488 : // METHODOLOGY EMPLOYED:
1489 : // A(i)*F(i,j)=A(j)*F(j,i); F(i,i)=0.; SUM(F(i,j)=1.0, j=1,N)
1490 : // Subroutine takes approximate view factors and enforces reciprocity by
1491 : // averaging AiFij and AjFji. Then it determines a set of row coefficients
1492 : // which can be multipled by each AF product to force the sum of AiFij for
1493 : // each row to equal Ai, and applies them. Completeness is checked, and if
1494 : // not satisfied, the AF averaging and row modifications are repeated until
1495 : // completeness is within a preselected small deviation from 1.0
1496 : // The routine also checks the number of surfaces and if N<=3, just enforces reciprocity.
1497 :
1498 : // Argument array dimensioning
1499 313 : EP_SIZE_CHECK(A, N);
1500 313 : F.dim(N, N);
1501 :
1502 : // SUBROUTINE PARAMETER DEFINITIONS:
1503 313 : Real64 constexpr PrimaryConvergence(0.001);
1504 313 : Real64 constexpr DifferenceConvergence(0.00001);
1505 :
1506 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
1507 : Real64 ConvrgNew;
1508 : Real64 CheckConvergeTolerance; // check value for actual warning
1509 :
1510 : bool Converged;
1511 : // OriginalCheckValue is the first pass at a completeness check. Even if this is zero,
1512 : // there is no guarantee that reciprocity is satisfied. As a result, the rest of this
1513 : // routine is needed to correct any issues even if the user defined view factors
1514 : // satisfy completeness (OriginalCheckValue = 0).
1515 313 : OriginalCheckValue = std::abs(sum(F) - N);
1516 :
1517 : // Allocate and zero arrays
1518 313 : Array2D<Real64> FixedAF(F); // store for largest area check
1519 :
1520 313 : Real64 ConvrgOld = 10.0;
1521 313 : Real64 LargestArea = maxval(A);
1522 313 : bool severeErrorPresent = false;
1523 : // Check for Strange Geometry
1524 : // When one surface has an area that exceeds the sum of all other surface areas in a zone,
1525 : // essentially the situation is a non-complete enclosure. As a result, the view factors
1526 : // calculated using the standard procedure below will not converge and may result in invalid
1527 : // view factors where either reciprocity or completeness is not satisfied. However, when
1528 : // the largest surface is just slightly smaller than the rest of the surface areas in the
1529 : // zone, it has been shown that there can still be problems. The correction below can
1530 : // be helpful in avoiding these problems. So, the criteria below (with the 0.99 term added
1531 : // into the comparison in the next line of code) intends to capture more cases that are
1532 : // "unbalanced" in surface area distribution so other strange cases can take advantage of
1533 : // this correction. The use of 0.99 is simply to provide some reasonable boundary numerically
1534 : // and does not have some derived theoretical basis.
1535 313 : if (LargestArea > 0.99 * (sum(A) - LargestArea) && (N > 3)) {
1536 1 : for (int i = 1; i <= N; ++i) {
1537 1 : if (LargestArea != A(i)) continue;
1538 1 : state.dataHeatBalIntRadExchg->LargestSurf = i;
1539 1 : break;
1540 : }
1541 1 : FixedAF(state.dataHeatBalIntRadExchg->LargestSurf, state.dataHeatBalIntRadExchg->LargestSurf) =
1542 1 : min(0.9, 1.2 * LargestArea / sum(A)); // Give self view to big surface
1543 : }
1544 :
1545 : // Set up AF matrix.
1546 313 : Array2D<Real64> AF(N, N); // = (AREA * DIRECT VIEW FACTOR) MATRIX
1547 2122 : for (int i = 1; i <= N; ++i) {
1548 13778 : for (int j = 1; j <= N; ++j) {
1549 11969 : AF(j, i) = FixedAF(j, i) * A(i);
1550 : }
1551 : }
1552 :
1553 : // Enforce reciprocity by averaging AiFij and AjFji
1554 313 : FixedAF = 0.5 * (AF + transpose(AF)); // Performance Slow way to average with transpose (heap use)
1555 :
1556 313 : AF.deallocate();
1557 :
1558 313 : Array2D<Real64> FixedF(N, N); // CORRECTED MATRIX OF VIEW FACTORS (N X N)
1559 :
1560 313 : NumIterations = 0;
1561 313 : RowSum = 0.0;
1562 : // Check for physically unreasonable enclosures.
1563 :
1564 313 : if (N <= 3) {
1565 195 : for (int i = 1; i <= N; ++i) {
1566 474 : for (int j = 1; j <= N; ++j) {
1567 337 : FixedF(j, i) = FixedAF(j, i) / A(i);
1568 : }
1569 : }
1570 :
1571 58 : ShowWarningError(state, format("Surfaces in Zone/Enclosure=\"{}\" do not define an enclosure.", enclName));
1572 116 : ShowContinueError(state, "Number of surfaces <= 3, view factors are set to force reciprocity but may not fulfill completeness.");
1573 116 : ShowContinueError(state, "Reciprocity means that radiant exchange between two surfaces will match and not lead to an energy loss.");
1574 116 : ShowContinueError(state,
1575 : "Completeness means that all of the view factors between a surface and the other surfaces in a zone add up to unity.");
1576 116 : ShowContinueError(state,
1577 : "So, when there are three or less surfaces in a zone, EnergyPlus will make sure there are no losses of energy but");
1578 116 : ShowContinueError(
1579 : state, "it will not exchange the full amount of radiation with the rest of the zone as it would if there was a completed enclosure.");
1580 :
1581 58 : RowSum = sum(FixedF);
1582 58 : if (RowSum > (N + 0.01)) {
1583 : // Reciprocity enforced but there is more radiation than possible somewhere since the sum of one of the rows
1584 : // is now greater than unity. This should not be allowed as it can cause issues with the heat balance.
1585 : // Correct this by finding the largest row summation and dividing all of the elements in the F matrix by
1586 : // this max summation. This will provide a cap on radiation so that no row has a sum greater than unity
1587 : // and will still maintain reciprocity.
1588 28 : Array1D<Real64> sumFixedF;
1589 : Real64 MaxFixedFRowSum;
1590 28 : sumFixedF.allocate(N);
1591 28 : sumFixedF = 0.0;
1592 96 : for (int i = 1; i <= N; ++i) {
1593 240 : for (int j = 1; j <= N; ++j) {
1594 172 : sumFixedF(i) += FixedF(i, j);
1595 : }
1596 68 : if (i == 1) {
1597 28 : MaxFixedFRowSum = sumFixedF(i);
1598 : } else {
1599 40 : if (sumFixedF(i) > MaxFixedFRowSum) MaxFixedFRowSum = sumFixedF(i);
1600 : }
1601 : }
1602 28 : sumFixedF.deallocate();
1603 28 : if (MaxFixedFRowSum < 1.0) {
1604 0 : ShowFatalError(state, " FixViewFactors: Three surface or less zone failing ViewFactorFix correction which should never happen.");
1605 : } else {
1606 28 : FixedF *= (1.0 / MaxFixedFRowSum);
1607 : }
1608 28 : RowSum = sum(FixedF); // needs to be recalculated
1609 28 : }
1610 58 : FinalCheckValue = FixedCheckValue = std::abs(RowSum - N);
1611 58 : F = FixedF;
1612 137 : for (int spaceNum : spaceNums) {
1613 79 : state.dataHeatBal->Zone(state.dataHeatBal->space(spaceNum).zoneNum).EnforcedReciprocity = true;
1614 : }
1615 58 : return; // Do not iterate, stop with reciprocity satisfied.
1616 :
1617 : } // N <= 3 Case
1618 :
1619 : // Regular fix cases
1620 255 : Array1D<Real64> RowCoefficient(N);
1621 255 : Converged = false;
1622 3369 : while (!Converged) {
1623 3114 : ++NumIterations;
1624 23950 : for (int i = 1; i <= N; ++i) {
1625 : // Determine row coefficients which will enforce closure.
1626 20836 : Real64 const sum_FixedAF_i(sum(FixedAF(_, i)));
1627 20836 : if (std::abs(sum_FixedAF_i) > 1.0e-10) {
1628 20836 : RowCoefficient(i) = A(i) / sum_FixedAF_i;
1629 : } else {
1630 0 : RowCoefficient(i) = 1.0;
1631 : }
1632 20836 : FixedAF(_, i) *= RowCoefficient(i);
1633 : }
1634 :
1635 : // Enforce reciprocity by averaging AiFij and AjFji
1636 3114 : FixedAF = 0.5 * (FixedAF + transpose(FixedAF));
1637 :
1638 : // Form FixedF matrix
1639 23950 : for (int i = 1; i <= N; ++i) {
1640 179108 : for (int j = 1; j <= N; ++j) {
1641 158272 : FixedF(j, i) = FixedAF(j, i) / A(i);
1642 158272 : if (std::abs(FixedF(j, i)) < 1.e-10) {
1643 29040 : FixedF(j, i) = 0.0;
1644 29040 : FixedAF(j, i) = 0.0;
1645 : }
1646 : }
1647 : }
1648 :
1649 3114 : ConvrgNew = std::abs(sum(FixedF) - N);
1650 3114 : if (std::abs(ConvrgOld - ConvrgNew) < DifferenceConvergence || ConvrgNew <= PrimaryConvergence) { // Change in sum of Fs must be small.
1651 255 : Converged = true;
1652 : }
1653 3114 : ConvrgOld = ConvrgNew;
1654 3114 : if (NumIterations > 400) { // If everything goes bad,enforce reciprocity and go home.
1655 : // Enforce reciprocity by averaging AiFij and AjFji
1656 0 : FixedAF = 0.5 * (FixedAF + transpose(FixedAF));
1657 :
1658 : // Form FixedF matrix
1659 0 : for (int i = 1; i <= N; ++i) {
1660 0 : for (int j = 1; j <= N; ++j) {
1661 0 : FixedF(j, i) = FixedAF(j, i) / A(i);
1662 : }
1663 : }
1664 0 : Real64 const sum_FixedF(sum(FixedF));
1665 0 : FinalCheckValue = FixedCheckValue = CheckConvergeTolerance = std::abs(sum_FixedF - N);
1666 0 : RowSum = sum_FixedF;
1667 0 : if (CheckConvergeTolerance > 0.005) {
1668 0 : if (CheckConvergeTolerance > 0.1) {
1669 0 : ShowSevereError(
1670 : state,
1671 0 : format("FixViewFactors: View factors convergence has failed and will lead to heat balance errors in zone=\"{}\".",
1672 : enclName));
1673 : }
1674 0 : ShowWarningError(
1675 : state,
1676 0 : format("FixViewFactors: View factors not complete. Check for bad surface descriptions or unenclosed zone=\"{}\".", enclName));
1677 0 : ShowContinueError(state,
1678 0 : format("Enforced reciprocity has tolerance (ideal is "
1679 : "0)=[{:.6R}], Row Sum (ideal is {})=[{:.2R}].",
1680 : CheckConvergeTolerance,
1681 : N,
1682 : RowSum));
1683 0 : ShowContinueError(state,
1684 : "If zone is unusual or tolerance is on the order of 0.001, view "
1685 : "factors might be OK but results should be checked carefully.");
1686 0 : if (anyIntMassInZone) {
1687 0 : ShowContinueError(state,
1688 : "For zones with internal mass like this one, this"
1689 : "can happen when the internal mass has an area that"
1690 : "is much larger than the other surfaces in the zone.");
1691 0 : ShowContinueError(state,
1692 : "If a single thermal mass element exists in this zone"
1693 : "that has an area that is larger than the sum of the"
1694 : "rest of the surface areas, consider breaking it up"
1695 : "into two or more separate internal mass elements.");
1696 : }
1697 : }
1698 0 : if (std::abs(FixedCheckValue) < std::abs(OriginalCheckValue)) {
1699 0 : F = FixedF;
1700 0 : FinalCheckValue = FixedCheckValue;
1701 : }
1702 0 : return;
1703 : }
1704 : }
1705 255 : FixedCheckValue = ConvrgNew;
1706 255 : if (FixedCheckValue < OriginalCheckValue) {
1707 0 : F = FixedF;
1708 0 : FinalCheckValue = FixedCheckValue;
1709 : } else {
1710 255 : FinalCheckValue = OriginalCheckValue;
1711 255 : RowSum = sum(FixedF);
1712 255 : if (std::abs(RowSum - N) < PrimaryConvergence) {
1713 255 : F = FixedF;
1714 255 : FinalCheckValue = FixedCheckValue;
1715 : } else {
1716 0 : ShowWarningError(
1717 : state,
1718 0 : format("FixViewFactors: View factors not complete. Check for bad surface descriptions or unenclosed zone=\"{}\".", enclName));
1719 : }
1720 : }
1721 255 : if (severeErrorPresent) {
1722 0 : ShowFatalError(state,
1723 : "FixViewFactors: View factor calculations significantly out "
1724 : "of tolerance. See above messages for more information.");
1725 : }
1726 429 : }
1727 :
1728 310 : bool DoesZoneHaveInternalMass(EnergyPlusData &state, int const numZoneSurfaces, const Array1D_int &surfPointer)
1729 : {
1730 2092 : for (int i = 1; i <= numZoneSurfaces; ++i) {
1731 1789 : if (state.dataSurface->Surface(surfPointer(i)).Class == DataSurfaces::SurfaceClass::IntMass) return true;
1732 : }
1733 303 : return false;
1734 : }
1735 :
1736 4062 : void CalcScriptF(EnergyPlusData &state,
1737 : int const N, // Number of surfaces
1738 : Array1D<Real64> const &A, // AREA VECTOR- ASSUMED,BE N ELEMENTS LONG
1739 : Array2<Real64> const &F, // DIRECT VIEW FACTOR MATRIX (N X N)
1740 : Array1D<Real64> &EMISS, // VECTOR OF SURFACE EMISSIVITIES
1741 : Array2<Real64> &ScriptF // MATRIX OF SCRIPT F FACTORS (N X N) //Tuned Transposed
1742 : )
1743 : {
1744 :
1745 : // SUBROUTINE INFORMATION:
1746 : // AUTHOR Curt Pedersen
1747 : // DATE WRITTEN 1980
1748 : // MODIFIED July 2000 (COP for the ASHRAE Loads Toolkit)
1749 : // RE-ENGINEERED September 2000 (RKS for EnergyPlus)
1750 : // RE-ENGINEERED June 2014 (Stuart Mentzer): Performance tuned
1751 :
1752 : // PURPOSE OF THIS SUBROUTINE:
1753 : // Determines Hottel's ScriptF coefficients which account for the total
1754 : // grey interchange between surfaces in an enclosure.
1755 :
1756 : // REFERENCES:
1757 : // Hottel, H. C. and A. F. Sarofim, Radiative Transfer, Ch 3, McGraw Hill, 1967.
1758 :
1759 : // SUBROUTINE ARGUMENTS:
1760 : // --Must satisfy reciprocity and completeness:
1761 : // A(i)*F(i,j)=A(j)*F(j,i); F(i,i)=0.; SUM(F(i,j)=1.0, j=1,N)
1762 :
1763 : // SUBROUTINE PARAMETER DEFINITIONS:
1764 4062 : Real64 constexpr MaxEmissLimit(0.99999); // Limit the emissivity internally/avoid a divide by zero error
1765 :
1766 : // Validate argument array dimensions
1767 4062 : assert(N >= 0); // Do we need to allow for N==0?
1768 4062 : assert((A.l() == 1) && (A.u() == N));
1769 4062 : assert((F.l1() == 1) && (F.u1() == N));
1770 4062 : assert((F.l2() == 1) && (F.u2() == N));
1771 4062 : assert((EMISS.l() == 1) && (EMISS.u() == N));
1772 4062 : assert(equal_dimensions(F, ScriptF));
1773 :
1774 : #ifdef EP_Count_Calls
1775 : ++state.dataTimingsData->NumCalcScriptF_Calls;
1776 : #endif
1777 :
1778 : // Load Cmatrix with AF (AREA * DIRECT VIEW FACTOR) matrix
1779 4062 : Array2D<Real64> Cmatrix(N, N); // = (AF - EMISS/REFLECTANCE) matrix (but plays other roles)
1780 4062 : assert(equal_dimensions(Cmatrix, F)); // For linear indexing
1781 4062 : Array2D<Real64>::size_type l(0u);
1782 31070 : for (int j = 1; j <= N; ++j) {
1783 214954 : for (int i = 1; i <= N; ++i, ++l) {
1784 187946 : Cmatrix[l] = A(i) * F[l]; // [ l ] == ( i, j )
1785 : }
1786 : }
1787 :
1788 : // Load Cmatrix with (AF - EMISS/REFLECTANCE) matrix
1789 4062 : Array1D<Real64> Excite(N); // Excitation vector = A*EMISS/REFLECTANCE
1790 4062 : l = 0u;
1791 31070 : for (int i = 1; i <= N; ++i, l += N + 1) {
1792 27008 : Real64 EMISS_i(EMISS(i));
1793 27008 : if (EMISS_i > MaxEmissLimit) { // Check/limit EMISS for this surface to avoid divide by zero below
1794 0 : EMISS_i = EMISS(i) = MaxEmissLimit;
1795 0 : ShowWarningError(state, "A thermal emissivity above 0.99999 was detected. This is not allowed. Value was reset to 0.99999");
1796 : }
1797 27008 : Real64 const EMISS_i_fac(A(i) / (1.0 - EMISS_i));
1798 27008 : Excite(i) = -EMISS_i * EMISS_i_fac; // Set up matrix columns for partial radiosity calculation
1799 27008 : Cmatrix[l] -= EMISS_i_fac; // Coefficient matrix for partial radiosity calculation // [ l ] == ( i, i )
1800 : }
1801 :
1802 4062 : Array2D<Real64> Cinverse(N, N); // Inverse of Cmatrix
1803 4062 : CalcMatrixInverse(Cmatrix, Cinverse); // SOLVE THE LINEAR SYSTEM
1804 4062 : Cmatrix.clear(); // Release memory ASAP
1805 :
1806 : // Scale Cinverse colums by excitation to get partial radiosity matrix
1807 4062 : l = 0u;
1808 31070 : for (int j = 1; j <= N; ++j) {
1809 27008 : Real64 const e_j(Excite(j));
1810 214954 : for (int i = 1; i <= N; ++i, ++l) {
1811 187946 : Cinverse[l] *= e_j; // [ l ] == ( i, j )
1812 : }
1813 : }
1814 4062 : Excite.clear(); // Release memory ASAP
1815 :
1816 : // Form Script F matrix transposed
1817 4062 : assert(equal_dimensions(Cinverse, ScriptF)); // For linear indexing
1818 4062 : Array2D<Real64>::size_type m(0u);
1819 31070 : for (int i = 1; i <= N; ++i) { // Inefficient order for cache but can reuse multiplier so faster choice depends on N
1820 27008 : Real64 const EMISS_i(EMISS(i));
1821 27008 : Real64 const EMISS_fac(EMISS_i / (1.0 - EMISS_i));
1822 27008 : l = static_cast<Array2D<Real64>::size_type>(i - 1);
1823 214954 : for (int j = 1; j <= N; ++j, l += N, ++m) {
1824 187946 : if (i == j) {
1825 : // ScriptF(I,J) = EMISS(I)/(1.0d0-EMISS(I))*(Jmatrix(I,J)-Delta*EMISS(I)), where Delta=1
1826 27008 : ScriptF[m] = EMISS_fac * (Cinverse[l] - EMISS_i); // [ l ] = ( i, j ), [ m ] == ( j, i )
1827 : } else {
1828 : // ScriptF(I,J) = EMISS(I)/(1.0d0-EMISS(I))*(Jmatrix(I,J)-Delta*EMISS(I)), where Delta=0
1829 160938 : ScriptF[m] = EMISS_fac * Cinverse[l]; // [ l ] == ( i, j ), [ m ] == ( j, i )
1830 : }
1831 : }
1832 : }
1833 4062 : }
1834 :
1835 4062 : void CalcMatrixInverse(Array2<Real64> &A, // Matrix: Gets reduced to L\U form
1836 : Array2<Real64> &I // Returned as inverse matrix
1837 : )
1838 : {
1839 : // SUBROUTINE INFORMATION:
1840 : // AUTHOR Jakob Asmundsson
1841 : // DATE WRITTEN January 1999
1842 : // MODIFIED September 2000 (RKS for EnergyPlus)
1843 : // RE-ENGINEERED June 2014 (Stuart Mentzer): Performance/memory tuning rewrite
1844 :
1845 : // PURPOSE OF THIS SUBROUTINE:
1846 : // To find the inverse of Matrix, using partial pivoting.
1847 :
1848 : // METHODOLOGY EMPLOYED:
1849 : // Inverse is found using partial pivoting and Gauss elimination
1850 :
1851 : // REFERENCES:
1852 : // Any Linear Algebra book
1853 :
1854 : // Validation
1855 4062 : assert(A.square());
1856 4062 : assert(A.I1() == A.I2());
1857 4062 : assert(equal_dimensions(A, I));
1858 :
1859 : // Initialization
1860 4062 : int const l(A.l1());
1861 4062 : int const u(A.u1());
1862 4062 : int const n(u - l + 1);
1863 4062 : I.to_identity(); // I starts out as identity
1864 :
1865 : // Could do row scaling here to improve condition and then check min pivot isn't too small
1866 :
1867 : // Compute in-place LU decomposition of [A|I] with row pivoting
1868 31070 : for (int i = l; i <= u; ++i) {
1869 :
1870 : // Find pivot row in column i below diagonal
1871 27008 : int iPiv = i;
1872 27008 : Real64 aPiv(std::abs(A(i, i)));
1873 27008 : int ik(A.index(i, i + 1));
1874 107477 : for (int k = i + 1; k <= u; ++k, ++ik) {
1875 80469 : Real64 const aAki(std::abs(A[ik])); // [ ik ] == ( i, k )
1876 80469 : if (aAki > aPiv) {
1877 0 : iPiv = k;
1878 0 : aPiv = aAki;
1879 : }
1880 : }
1881 27008 : assert(aPiv != 0.0); //? Is zero pivot possible for some user inputs? If so if test/handler needed
1882 :
1883 : // Swap row i with pivot row
1884 27008 : if (iPiv != i) {
1885 0 : int ji(A.index(l, i)); // [ ji ] == ( j, i )
1886 0 : int pj(A.index(l, iPiv)); // [ pj ] == ( j, iPiv )
1887 0 : for (int j = l; j <= u; ++j, ji += n, pj += n) {
1888 0 : Real64 const Aij(A[ji]);
1889 0 : A[ji] = A[pj];
1890 0 : A[pj] = Aij;
1891 0 : Real64 const Iij(I[ji]);
1892 0 : I[ji] = I[pj];
1893 0 : I[pj] = Iij;
1894 : }
1895 : }
1896 :
1897 : // Put multipliers in column i and reduce block below A(i,i)
1898 27008 : Real64 const Aii_inv(1.0 / A(i, i));
1899 107477 : for (int k = i + 1; k <= u; ++k) {
1900 80469 : Real64 const multiplier(A(i, k) * Aii_inv);
1901 80469 : A(i, k) = multiplier;
1902 80469 : if (multiplier != 0.0) {
1903 77196 : int ji(A.index(i + 1, i)); // [ ji ] == ( j, i )
1904 77196 : int jk(A.index(i + 1, k)); // [ jk ] == ( j, k )
1905 427655 : for (int j = i + 1; j <= u; ++j, ji += n, jk += n) {
1906 350459 : A[jk] -= multiplier * A[ji];
1907 : }
1908 77196 : ji = A.index(l, i);
1909 77196 : jk = A.index(l, k);
1910 649752 : for (int j = l; j <= u; ++j, ji += n, jk += n) {
1911 572556 : Real64 const Iij(I[ji]);
1912 572556 : if (Iij != 0.0) {
1913 222097 : I[jk] -= multiplier * Iij;
1914 : }
1915 : }
1916 : }
1917 : }
1918 : }
1919 :
1920 : // Perform back-substitution on [U|I] to put inverse in I
1921 31070 : for (int k = u; k >= l; --k) {
1922 27008 : Real64 const Akk_inv(1.0 / A(k, k));
1923 27008 : int jk(A.index(l, k)); // [ jk ] == ( j, k )
1924 214954 : for (int j = l; j <= u; ++j, jk += n) {
1925 187946 : I[jk] *= Akk_inv;
1926 : }
1927 27008 : int ik(A.index(k, l)); // [ ik ] == ( i, k )
1928 107477 : for (int i = l; i < k; ++i, ++ik) { // Eliminate kth column entries from I in rows above k
1929 80469 : Real64 const Aik(A[ik]);
1930 80469 : int ji(A.index(l, i)); // [ ji ] == ( j, i )
1931 80469 : int jm(A.index(l, k)); // [ jm ] == ( k, j )
1932 676110 : for (int j = l; j <= u; ++j, ji += n, jm += n) {
1933 595641 : I[ji] -= Aik * I[jm];
1934 : }
1935 : }
1936 : }
1937 4062 : }
1938 :
1939 3 : void CalcFMRT(EnergyPlusData &state,
1940 : int const N, // Number of surfaces
1941 : Array1D<Real64> const &A, // AREA VECTOR- ASSUMED,BE N ELEMENTS LONG
1942 : Array1D<Real64> &FMRT // VECTOR OF MEAN RADIANT TEMPERATURE "VIEW FACTORS"
1943 : )
1944 : {
1945 3 : double sumAF = 0.0;
1946 10 : for (int iS = 0; iS < N; iS++) {
1947 7 : FMRT[iS] = 1.0;
1948 7 : sumAF += A[iS];
1949 : }
1950 :
1951 3 : constexpr int maxIt = 100;
1952 3 : constexpr double tol = 0.0001;
1953 : double fLast;
1954 3 : double sumAFNew = sumAF;
1955 11 : for (unsigned i = 0; i < maxIt; i++) {
1956 11 : double fChange = 0.;
1957 11 : bool errorsFound(false);
1958 11 : sumAF = sumAFNew;
1959 11 : sumAFNew = 0.0;
1960 33 : for (int iS = 0; iS < N; iS++) {
1961 23 : fLast = FMRT[iS];
1962 23 : FMRT[iS] = 1. / (1. - A[iS] * FMRT[iS] / (sumAF));
1963 23 : if (FMRT[iS] > 100.) {
1964 1 : errorsFound = true;
1965 2 : ShowSevereError(state, "Geometry not compatible with Carroll MRT Zone Radiant Exchange method.");
1966 1 : break;
1967 : }
1968 22 : fChange += fabs(FMRT[iS] - fLast);
1969 22 : sumAFNew += A[iS] * FMRT[iS];
1970 : }
1971 :
1972 11 : if (errorsFound || fChange / N < tol) {
1973 : break;
1974 : }
1975 8 : if (i >= maxIt) {
1976 0 : errorsFound = true;
1977 0 : ShowSevereError(state, "Carroll MRT Zone Radiant Exchange method unable to converge on \"view factor\" calculation.");
1978 : }
1979 8 : if (errorsFound) {
1980 0 : ShowFatalError(state, "CalcFMRT: Errors found while calculating mean radiant temperature view factors. Program terminated.");
1981 : }
1982 : }
1983 3 : return;
1984 : }
1985 :
1986 3 : void CalcFp(int const N, // Number of surfaces
1987 : Array1D<Real64> const &EMISS, // VECTOR OF SURFACE EMISSIVITIES
1988 : Array1D<Real64> const &FMRT, // VECTOR OF MEAN RADIANT TEMPERATURE "VIEW FACTORS"
1989 : Array1D<Real64> &Fp // VECTOR OF OPPENHEIM RESISTANCE VALUES
1990 : )
1991 : {
1992 10 : for (int iS = 0; iS < N; iS++) {
1993 7 : Fp[iS] = Constant::StefanBoltzmann * EMISS[iS] / (EMISS[iS] / FMRT[iS] + 1. - EMISS[iS]); // actually sigma *
1994 : }
1995 3 : return;
1996 : }
1997 :
1998 11 : int GetRadiantSystemSurface(EnergyPlusData &state,
1999 : std::string const &cCurrentModuleObject, // Calling Object type
2000 : std::string const &RadSysName, // Calling Object name
2001 : int const RadSysZoneNum, // Radiant system zone number
2002 : std::string const &SurfaceName, // Referenced surface name
2003 : bool &ErrorsFound // True when errors are found
2004 : )
2005 : {
2006 : static constexpr std::string_view routineName("GetRadiantSystemSurface: "); // include trailing blank space
2007 :
2008 : // For radiant zone equipment, find the referenced surface and check if it is in the same zone or radiant enclosure
2009 11 : int const surfNum = Util::FindItemInList(SurfaceName, state.dataSurface->Surface);
2010 :
2011 : // Trap for surfaces that do not exist
2012 11 : if (surfNum == 0) {
2013 0 : ShowSevereError(state, format("{}Invalid Surface name = {}", routineName, SurfaceName));
2014 0 : ShowContinueError(state, format("Occurs for {} = {}", cCurrentModuleObject, RadSysName));
2015 0 : ErrorsFound = true;
2016 0 : return surfNum;
2017 : }
2018 :
2019 11 : if (RadSysZoneNum == 0) {
2020 0 : ShowSevereError(state, format("{}Invalid Zone number passed by {} = {}", routineName, cCurrentModuleObject, RadSysName));
2021 0 : ErrorsFound = true;
2022 0 : return surfNum;
2023 : }
2024 :
2025 : // Check if the surface and equipment are in the same zone
2026 11 : int const surfZoneNum = state.dataSurface->Surface(surfNum).Zone;
2027 11 : if (surfZoneNum == 0) {
2028 : // This should never happen
2029 0 : ShowSevereError(state,
2030 0 : format("{}Somehow the surface zone number is zero for{} = {} and Surface = {}",
2031 : routineName,
2032 : cCurrentModuleObject,
2033 : RadSysName,
2034 : SurfaceName)); // LCOV_EXCL_LINE
2035 : ErrorsFound = true; // LCOV_EXCL_LINE
2036 11 : } else if (surfZoneNum != RadSysZoneNum) {
2037 0 : ShowSevereError(state, format("{}Surface = {} is not in the same zone as the radiant equipment.", routineName, SurfaceName));
2038 0 : ShowContinueError(state, format("Surface zone or enclosure = {}", state.dataHeatBal->Zone(surfZoneNum).Name));
2039 0 : ShowContinueError(state, format("Radiant equipment zone or enclosure = {}", state.dataHeatBal->Zone(RadSysZoneNum).Name));
2040 0 : ShowContinueError(state, format("Occurs for {} = {}", cCurrentModuleObject, RadSysName));
2041 0 : ErrorsFound = true;
2042 : }
2043 11 : return surfNum;
2044 : }
2045 :
2046 : } // namespace HeatBalanceIntRadExchange
2047 :
2048 : } // namespace EnergyPlus
|
