### Abstract

Research shows that real-time visual feedback can greatly assist the process of identifying and understanding complex cause and effect relationships. Whilst not all physical processes in the performance analysis of buildings can be simulated in real-time, some can. By investigating these and combining them with modern games technology, it is possible to create highly visual and interactive design tools that allow all the governing parameters of a simulation, including geometry, to be manipulated by the user in real-time with detailed contextual results updating dynamically with each interactive change.

This paper introduces a web-based daylighting simulation tool that does just that, using variants of the daylight coefficients and split-flux methods implemented on the GPU to calculate the real-time spatial distribution of daylight factors across a simple rectangular room. This tool\’s purpose is primarily educational, allowing users to gain a comprehensive understanding of the relationships between room dimensions, window placement and daylight distribution through a process of deliberate investigative play. A detailed parametric comparison with Radiance simulations resulted in a small modifier to the implemented method that produces a robust correlation which the authors argue makes this tool a valuable educational resource with potential applications in early design decision-making.

# Introduction

This work is part of ongoing research into the development of design tools that utilise analytical calculations performed on the graphics processing unit (GPU) to create highly dynamic and interactive building simulation environments. Thus the primary goal is to undertake analysis that is fast enough to provide dynamic visual feedback during model manipulations and accurate enough for that feedback to be meaningful and useful.

This paper describes a prototype web-based tool that utilises WebGL and custom shaders to dynamically calculate the spatial distribution of daylight within a rectangular room and update it in real-time as the user interactively manipulates room dimensions, surface properties, window sizes and their positions within the envelope.

The choice of daylight calculation method and its implementation within the prototype is described in detail as well as the process of comparing results against those of the same room configurations simulated in Radiance. Whilst validation is important to show that the approach taken does not yield misleading results, the main aim of the comparative analysis was to gain a better understanding of where the implemented method diverges from others and why.

The intention is for that information, and the experience and insights gained through the initial implementation, to inform future work by the authors and others in the development of much faster and even more accurate daylight algorithms that are better able to exploit the capabilities of modern GPUs. It will also serve to expose those areas of the initial method that could be improved or approached differently to maximise convergence and accommodate more sophisticated requirements such as complex glazing systems, dynamic shading devices and time-based variation in sky conditions.

# The Spatial Model

The geometric model used in this prototype tool is limited to a simple rectangular room with no internal obstructions, but with adjustable wall thickness and any number of rectangular windows in any wall. The width and depth of window frames are also customisable and can include any number of mullions and transoms. Some example screenshots of the room model and its controls are shown in Figure 1.

The decision to use such a simple model in the initial implementation was made for the following reasons:

• The simplified method used for calculating the internally reflected component is best suited to convex spaces with no internal obstructions (Winkelmann and Selkowitz, 1985),

• An axially-aligned rectangular plan shape greatly simplifies both the development and potential intuitiveness of interactive user manipulation techniques for the room and its windows,

• A rectangular plan shape also simplifies the process of automatically fitting an analysis grid over the work-plane and within the internal bounds of the space,

• Whilst not representative of all possible room shapes and sizes, such a configuration is broadly representative of a wide range of actual design conditions likely to be faced by designers,

• The potential ease with which the room can be manipulated to match a given design condition means that the tool does not require an extensive model creation or geometry import phase before data can be generated and displayed, and

• The simplicity of the geometry allows for experimentation with implicit design rules such as not allowing windows to overlap whilst dragging and how to adjust windows when walls are interactively reduced beyond the required size to accommodate them.

It is intended that future tools will better accommodate the detailed geometry and specific materiality of actual design spaces and allow CAD/BIM data to be imported. However, as a proof-of-concept, this tool is a necessary first step towards that goal.

# Available Daylight Methods

There are basically three methods for simulating the spatial distribution of diffuse illuminance within a room; these being the split-flux, radiosity and ray-tracing methods. The simplest of these is the split-flux method (Hopkinson et al, 1954), (Lynes, 1968) which splits incoming illuminance into downward and upward flux components, distributing them evenly over appropriate internal surfaces and applying area-weighted reflectance and absorptance. After this first distribution, average flux balance is used to determine statistically the total flux lost to surface absorption or back outside through a window.

The radiosity method (Goral et al, 1984) uses radiant energy transfer and form factors to determine visibility between surfaces and therefore both surface illuminance and luminance. The direct energy from all light sources is first distributed throughout the space and then radiant exchange between surfaces is calculated iteratively until all is absorbed. This gives the luminous distribution of each surface from which the illuminance at any point can calculated.

The ray-tracing method (Whitted, 1979) determines the progress of rays through a space, as they hit and reflect off surfaces within it. Diffuse luminance/illuminance is calculated for each visible point within an image by distributing rays back into the space to see which light sources and other illuminated surfaces are \‘visible\’ to that point.

Of these, the split-flux method is the simplest to implement and the fastest to calculate. The radiosity method is generally faster than ray-tracing, but often only when multiple views are being generated as the radiant exchange only needs to be solved once whereas extensive ray-tracing is required for each view. However, both radiosity and ray-tracing are resource intensive and orders of magnitude slower than the split-flux method.

One of the core design decisions in the development of this tool was to start by adapting a widely used and well understood methodology so that it could be easily tested and validated throughout each step of the implementation. As the split-flux method is widely used for daylight estimation at early design stage, and is the default daylight calculation method used within EnergyPlus (DOE, 2017), it was selected for initial implementation within this work.

# Daylight Calculation Method

The very first implementation iteration was actually based on the BRE Daylight Factor Protractor method (BRE, 1986). This method was particularly interesting because the values in each protractor already account for the effects of angular-dependent transmission through glass, the luminous distribution of an overcast sky and the conversion of aperture solid angle to daylight factor contribution. Thus, by digitising the daylight factor protractor and converting that data to high-resolution look-up tables, a significant number of quite complex calculations can be sidestepped.

The daylight factor at a point inside a room is basically the ratio of total received illuminance to the total available illuminance from an unobstructed overcast sky, expressed as a percentage. The total received illuminance within the room model is the sum of three components:

• Sky Component (SC):
Light passing through a window directly from the sky,

• Externally Reflected Component (ERC):
Light passing through a window after reflecting off outside surfaces or external obstructions, and

• Internally Reflected Component (IRC):
Light reflected from internal surfaces after bouncing around inside the space.

The daylight factor protractor method can be used to calculate both the sky and externally reflected components, whilst the split-flux method is used to determine the internally reflected component.

As described in the BRE method, the following equation is used for the room-averaged internally reflected component:

$$IRC_{avg} = \frac{0.85 . W}{A(1-R)} \times (C.R_{fw} + 5.R_{cw}) \tag{1}$$

Where:

• W = The total window area in (m^2^),

• A = The total area of internal surfaces in the room (m^2^).

• R = The area-weighted reflectance of internal surfaces in the room (0-1),

• C = A coefficient based on the average altitude to the top of external obstructions as seen and measured from the window center, read from a table.

• Rfw = The area-weighted reflectance of the floor and that area of adjacent and opposite walls below the center of the window (not including the window wall).

• Rcw = The area-weighted reflectance of the ceiling and that area of adjacent and opposite walls above the center of the window (not including the window wall).

• 0.85 = The average diffuse transmittance of a single pane of clear float glass, and

• 5 = The default reflectance of the ground plane and external obstructions, given as a percentage.

The room model used in this tool allows windows to be placed at any height within a wall and even stacked one above the other, which means that the variables C, Rfw and Rcw could vary quite significantly from window to window even in the same wall. Thus, rather than calculate some kind of average center height for windows across the whole room, this implementation assumes that diffuse daylight entering each window is incoherent and is therefore additive. Equation 1 is then solved once for each window, using window-specific values for those variables and substituting each window\’s surface area for W. The IRC contribution of each window is then summed to give the total IRC for the room.

## Using the GPU

The look-up tables for daylight factor are easily encoded into two dimensional textures. The GPU was then used to calculate the plan and sectional angles of each aperture from each point on the spatial daylight grid. These angles were used to linearly interpolate between looked-up texture values to derive the sky component of the daylight factor at each point.

Straightforward trigonometric calculations such as these are trivial for most modern GPUs and they are able to process grid points in parallel. Quantifying the exact capacity for parallelisation on any GPU is a fraught process as there are many interdependent criteria as well as different chip architectures and therefore different nomenclatures. However as a rough comparative measure, nVidia GeForce GPU chips range from around 700 parallel units/cores to as many as 4000 at the high end. AMD Radeon GPUs range from around 1,000 to nearly 2,500 whilst the Intel HD 500 series range from 110 to 1100 (TechPowerUp, 2017). Even the GPUs found in most phones and tablets are able to process more than a hundred points at once. At 250mm centers, a 6x10m room would contain around 960 grid points over the work-plane.

Profiling the execution of the tool during interactive user manipulation of a 6x10m room shows that a full recalculation and update of the model and user interface takes around 12-15 milliseconds on a standard i7 2014 MacBook Pro, and around 28-36 milliseconds on an iPad Air and Galaxy Note 4 phone. The actual daylight calculations across the grid take up less than 3% of that time on the MacBook Pro and only 4% on the iPad and Galaxy Note. More importantly, the majority of that 3 to 4% is actually spent encoding the updated room and aperture data into shader uniforms on the CPU and then decoding the GPU results back to grid array values. The remaining 96 to 97% of that CPU and GPU time is spent regenerating and re-tessellating the polygons forming room geometry and daylight grid contours, as well as updating the various graphical user interface elements.

Thus, this work has found that this approach to daylighting in such a simple model can be so highly optimised that it becomes an insignificant component of the frame-by-frame workload required to visualise changes in real-time.

Whilst this was something of a surprise, even more surprising was the fact that the same approach can be coded entirely in Javascript, without using the GPU at all, and the execution times are very similar. The dynamic code optimisation capabilities of Javascript compilers in most modern browsers means that highly repetitive numeric calculations are quickly optimised to near native speeds. Thus the time spent by the CPU decoding the GPU results is roughly equivalent to it actually performing all the calculations itself.

All of this essentially meant that there was significant capacity for additional computation, on both the GPU or CPU, before the frame rates of dynamic updates were noticeably impacted. How much capacity is obviously device-dependant, but even the iPad and Galaxy Note were able to sustain frame rates of more than 20fps during an interactive drag.

## Switching to Daylight Coefficients

This additional computational capacity allowed some of the limitations of the protractor-based method to be reconsidered. As protractors have only been published for a very specific set of sky luminance distributions and glazing types, handling the full 16 CIE Standard General Sky types for climate based daylight modelling (CBDM) and accommodating more complex glazing and shading systems requires a different approach.

Recent work on CBDM has used hemispheric subdivision techniques to simulate skylight and sunlight as a series of individual sky patches, each with varying luminance (Bourgeois D, Reinhart CF, Ward G. 2008). Replacing the protractor-based look-up tables with an array of sky patches and then using the GPU to determine which patches are visible from each grid point has several advantages:

• The horizontal and vertical angles of each sky patch can be easily encoded as a two dimensional texture and sent to the GPU in the same way as look-up tables,

• Patch visibility can be cached for each grid point, along with additional information such as the aperture through which it is visible or a reference to any internal or external obstructing surface(s),

• The sky dome is considered to be sufficiently distant that the angles of each patch are the same for each grid point, meaning that angles of transmission through glazing or surface incidence do not need to be cached as they can be easily determined by look-up as required,

• Knowing surface intersection and/or aperture transmission angles for each patch allows more complex glazing, shading and light redirection systems to be analysed and even BRDF/BSDF functions incorporated,

• Cached aperture and obstruction data can be used to further optimise interactive manipulations - such as changing window frame configurations or resizing an aperture which requires only those patches previously passing through that aperture or obstructed by the containing wall to recalculated,

• Having each grid point check the same number of patches and therefore reference the same texture coordinates at the same time and in the same order is a process particularly suitable for the single instruction, multiple data (SIMD) architectures of almost all modern GPUs, making it fast and efficient even when several thousand sky patches are used, and

• Unlike the protractor-based approach, the sky subdivision method and number of sky patches can be dynamically tailored to both the complexity of the model and the capabilities of the hardware on which it is running, even during an interactive event if updates are found to be taking too long. This means potentially switching between equal-angle, equal-area or the Tregenza/Reinhart distribution, with different subdivision angles.

Figure 3 shows a simple example of how sky patches work. One of the sky subdivision methods is used to divide the sky into patches. The visibility of each sky patch is then calculated and stored at each grid point on the work plane, along with metadata such as the aperture passed through or the internal/external surface obstruction. The instantaneous or cumulative sky luminance distribution is then mapped over the patches and the contributions of all visible or semi-visible patches at each point are summed to determine the direct illuminance. In actuality it is slightly more complex in that the sky patch visibility can also include surface and glazing incidence effects and the sky luminance distribution is often split into direct and diffuse components and considered separately.

# Comparison of Results

Both the BRE and EnergyPlus versions of the split-flux method are all based on early work by Lynes (1979) on the derivation of a daylight factor formula for side-lit rectangular rooms. Some researchers have found that Lynes\’ formula and its close derivatives have a useful correlation with both measured daylight data (Crisp, Littlefair, 1984) and Radiance simulations of the same spaces (Reinhart, 2010). Detailed studies by Versage et al. (2010) and Yoon et al. (2014) indicate that, for south oriented windows, the split-flux method predicts higher illuminances than the radiosity and ray-tracing algorithms as the distance from windows increases. However, follow-up studies by Yoon et al. (2014) comparing other orientations indicate that this is not always the case and that much depends on the settings used in each radiosity or ray-tracing run.

All of these studies comparing the split-flux with other methods have focused on absolute accuracy and differences in individual values. However, the aim of this tool is to identity and illustrate relationships, so relative accuracy when calculation parameters change is of more importance. For example, does predicted daylight fall by the same percentage in all methods when the window area is halved? If the relative trends match and overall correlation is high, then the results can still provide meaningful and useful insight and design guidance even if the absolute values do not exactly match.

Radiance (Ward, and Rubinstein, 1988) is a widely used and highly validated daylight simulation program (Mardeljevic, 1997) developed by Greg Ward and Lawrence Berkeley National Laboratories. It is based on a variant of the ray-tracing method and is used at all levels of lighting and daylighting design as the reference simulation tool.

To investigate both the absolute and relative accuracy, as well as the overall correlation, results from the GPU-based split-flux method were compared with those from spaces with exactly the same configuration simulated in Radiance. As comprehensive datasets of measured light levels in real rooms are limited (Osborne and Donn, 2011), and those datasets mainly exist as a result of having been used to validate Radiance (Mardaljevic, 2000), validating this tool against Radiance allows for comparison over a much wider range of room sizes and aperture layouts than would be possible using only measured data.

Also, in order to ensure that any correlation or otherwise was not simply a matter of configuration coincidence, a parametric comparison was undertaken over a range of room and aperture sizes, aperture positions, frame sizes, surface properties and work-plane heights. This parametric comparison looked not only at absolute differences in value, but at relative trends in the way the two daylight distributions changed as each parameter was varied.

## Room Configurations and Parametrics

The default room settings used for all runs in which those parameters were not being varied were as follows:

• An internal room size of 5000 mm in the X axis, 7500 mm on the Y axis, a height of 2400mm in the Z axis and a wall thickness of 210mm,

• Aperture sizes with a height of 1200 mm and widths of 2000 mm in the short wall and 3000mm in the longer wall,

• Window frames 50 mm wide and 100 mm deep, centered within the 210 mm wall,

• Surface reflectances of 0.6 for walls and window frames, 0.7 for the ceiling and 0.4 for the floor, and

• A horizontal work-plane placed 750 mm above the floor level.

Based on these defaults, five basic aperture configurations were used:

• A single window in the shorter wall,

• A single window in the longer wall,

• Two adjacent windows, one in the shorter and one in the longer wall,

• Three windows, one in each of the two shorter walls and one in the longer wall, and

• Four windows, one in each wall.

For each aperture configuration, variations in the following parameters were then studied:

• Aperture Width: Varied from 250 mm wide to the full width of the containing wall in increments of 250 mm.

• Aperture Height: Heights were varied from 250 mm to full wall width in increments of 250 mm,

• Aperture Sill Height: Varied from 50mm to 1350 mm in increments of 250 mm,

• Frame Width: Varied from 0 (no frame) to 500 mm in increments of 50 mm.

• Surface Reflectance: Varied from 0.2 to 0.8 for walls and window frames, 0.3 to 0.9 for ceilings and 0.1 to 0.7 for the floor.

• Work-Plane Height: Varied from 250 mm above the floor to 2250 mm in increments of 250 mm.

• Room Size: The width and depth were varied separately from 2500 mm to 20000 mm in increments of 500 mm.

Generation of the parametric models and automation of the Radiance runs for the validation was performed using a Node.js (Node.js Foundation, 2017) script. As both the browser app itself and the node script were written in Javascript, they each share the same JSON room description and much of the same geometry code. The difference being that the browser app generates indexed triangle and line buffers required by WebGL whilst the node script generates polygonal geometry required for Radiance input files.

The basic steps involved in generating a Radiance run with a given room configuration and an example base name of PROJECT are as follows:

1. Generate materials and model geometry - create PROJECT.rad file,

2. Generate ground definition and an overcast sky description with an overall horizontal irradiance value of 100 W/m2 - create PROJECT\_sky.rad file,

3. Generate daylight grid over the workplane, saving positions and upwards-facing normal vectors - create PROJECT.pts file,

4. Generate run settings and optional sectional plan view - create PROJECT.rif file,

5. Invoke rad with the PROJECT.rif to create the octree, ambient and optional sectional plan view image files,

6. Pipe daylight grid points through rtrace to calculate horizontal irradiance values - create PROJECT.dat file,

7. Read completed PROJECT.dat file, converting irradiance back to illuminance and marrying the resulting values with their appropriate daylight grid points, and

8. Compare Radiance values over the work-plane with those from the browser app to calculate the correlation coefficient, root mean squared error, average difference and standard deviation.

#### Step 1: Generating Radiance Materials

The internal floor, wall and ceilings within the room model are each defined by their surface reflectance and the Lynes daylight factor formula assumes that all internal reflectance is Lambertian. To model this in Radiance, a plastic material type was used which requires parameters for its RGB surface reflectance, specularity and surface roughness. As the desired result of this analysis is daylight factor, a gray color for all surfaces is used with a low specularity (0.05) and medium surface roughness (0.1) to ensure a matte finish.

Thus, for a wall with a surface reflectance of 0.6, the Radiance material output would be:

    void plastic wall
0
0
5 0.6 0.6 0.6 0.05 0.1


For glazing, a Radiance glass material was used. This type of material is optimized for thin glass surfaces without internal reflections and the required parameters are its RGB transmissivity at normal incidence (tn). Within the room model, glazing is defined by its transmittance (Tn) which is the value published by most manufacturers, which the Radiance reference manual (Berkeley Lab, 2017) recommends converting to transmissivity by the following equation:

$$t_n = \frac{\sqrt{(0.8402528435 + (0.0072522239 . T_n^2)) - 0.9166530661}} {0.0036261119 / T_n} \tag{2}$$

Thus, for glazing with a transmittance of 0.6, the Radiance material output would be:

    void glass glazing
0
0
3 0.65405 0.65405 0.65405


#### Step 2: Ground and Overcast Sky Conditions

To define an overcast sky in Radiance, the gensky command is used. To ensure that the Sun is not affecting the zenith brightness, even in a cloudy sky with no direct sunlight, a solar altitude of 60 degrees and an azimuth of 0 are explicitly set using the -ang option. This means that no date, time or geographic location settings are required. The cloudy sky is specified with the -c option.

The -B option is also specified with a value of 100.0, which is actually quite important as this explicitly sets the horizontal diffuse irradiance in W/m^2^ that the generated sky must provide. As we are comparing daylight factors, setting the external unobstructed horizontal diffuse irradiance to 100 W/m^2^ means that the calculated internal horizontal diffuse irradiance values at each grid point will automatically be daylight factors given as percentages, without requiring any conversion. The Radiance command to generate the required sky is then simply:

    !gensky -ang 90.0 0.0 -c -B 100.0


This command generates a skyfunc modifier which must then be applied to the model. This in done in Radiance by first defining a glow material that passes through the unmodified values from skyfunc hence the first three RGB values of 1. The fourth parameter is the maximum radius for shadow testing which, as this is not required for diffuse skylight, is set to zero and ignored.

    skyfunc glow sky_mat
0
0
4 1 1 1 0


This glow is then assigned to a directional light source positioned directly above the model. The first three parameters give the direction vector towards the source. The fourth is the angle of light distribution from the source which, in this case, is a full 180 degrees.

    sky_mat source sky
0
0
4 0 0 1 180


Defining a ground surface is done in the same way as the sky, but by modifying values from skyfunc by the required ground reflectance and positioning the source below the model using a direction vector that points downwards. Whilst the internally reflected component calculated with Equation 1 assumes a ground reflectance of 5%, when a value of 0.05 is used in Radiance it all but eliminates ground effects. Thus, in this comparison work a value of 0.2 is used for the first three parameters of the ground glow modifier to more realistically simulate average ground conditions.

    skyfunc glow ground_glow
0
0
4 0.2 0.2 0.2 0

ground_glow source ground
0
0
4 0 0 -1 180


Instead of defining a ground source, in theory it is possible to use the -g option within the gensky command to set the ground reflectance. In this case, any upward sky vectors processed by skyfunc are modified by this value so no ground source is required. However, this produces results very close to those with no ground reflectance at all and which do not respond to changes in ground reflectance value. The ground source approach described above does respond to changes in ground reflectance value and increases internal illuminance values as expected.

#### Step 3: Generating Grid Points

To ensure that grid points in each calculation are taken at the same positions, the JSON room description generated by the browser app includes a grid.points array containing the absolute X, Y and Z position of each grid point on the work-plane. If this property is not present in the room description, grid point positions are generated using the same algorithm (and the same shared code) as the browser app.

It is important to note that each grid point is located at the center of a rectangular region of the work-plane and represents the average value over that region. Thus, even though points at the edge of the grid appear spaced off the walls by half an increment, they do represent the entire work-plane area right up to each wall as their region extends half an increment in each direction.

A design decision was made to visually represent the daylight grid as smoothly interpolated triangles drawn between each grid point, rather than as a series of non-interpolated coloured rectangles for each representative region. Whilst non-interpolated coloured regions better represent the underlying spatial sampling used within the tool, smoothly interpolated triangles better represent the actual nature of the data being sampled.

Grid points are added to the .pts file one per line with space separated X, Y and Z values followed by space separated NX, NY and NZ values defining the unit surface normal of the work-plane. In this comparison, the work-plane was always horizontal so the normal values were 0 0 1.

#### Step 4: Generating Run Settings

One important lesson learnt as part of this work was that the ambient file is governed entirely by the settings used to create it. Thus, it does not appear possible to create an initial ambient file with low resolution ambient settings and then generate additional renders with higher resolution settings and expect improved results. The sampling and distribution of points in the ambient file appear to stay the same as when it was first created.

This is significant as the first approach used in this work was to use an OPTFILE entry in the .rif file generated in this step, which causes Radiance to create a .opt file with all the detailed rpict options it chose during the first run in Step 5. Then, as part of Step 6, this .opt file was read in and the parameter values compared against a range of model metrics to determine if any needed tweaking before the grid data point values were then calculated.

As the grid point data is only ever viewed within the web-based app and never directly seen within the context of the Radiance rendered model, it took some time to realise that the results were not responsive to changes in certain parameters, in this case the ambient parameters -aa, -ab, -ad, -ar and -as. This also acted to mask much of the responsiveness to other parameters.

Thus, the only approach that works is to analyse room model metrics during this step, prior to any Radiance runs, and determine the appropriate settings for each important rpict parameter before adding them to the “render = “ line of the .rif file. This allows Radiance to use the same parameter set for rpict in Step 5 and rtrace in Step 6.

However, this requires not only a comprehensive understanding of all the rpict/rtrace parameters and what they do, but significant for-knowledge of how they are likely to affect the model simulation and interact with other parameters. This approach was viable in this particular case as it is based on a relatively simple room model which was subject to extensive parametric testing. However the authors suspect that it will become less viable when extended to arbitrary room geometry and material properties in future versions.

##### Criteria for Changing Settings

The fundamental criteria for determining if parameter values need tweaking is the appearance of splotches within the rendered Radiance image. This indicates that the current ambient settings are not sufficient to deal with the gradients of change occurring within the lighting environment and insufficient rays are being generated to detect all the light sources in the same way across each surface. If the parameters are not tweaked, this will manifest itself as small apparently random variations in grid point results, making the work-plane daylight distribution appear slightly bumpy and reducing correlation with the smoother results from the simpler method.

As shown in Figure 4, the occurrence of splotches increases significantly as the aperture area becomes quite small compared to the total internal surface area of the room. This can occur in regularly sized rooms with very small windows or in very large rooms with a comparatively low ceiling height. This work has found that this becomes noticeable when total aperture area is below around 7.5% of the total internal surface area.

To solve this problem, many Radiance users switch to modelling each aperture as an illum source (Berkeley Lab, 2017), in which case each aperture is treated as an artificial light and the mkillum utility is used to pre-calculate their illuminances from the current sky conditions.

The alternative is to significantly increase ambient settings to compensate, which results in a commensurately significant increase in calculation time. When using the highest settings shown on the right in Figure 5, calculation times increase from around 20 seconds on a i7 2014 MacBook Pro to around three (3) hours.

##### Rendering a Redundant Image

If no view is specified, the rad command will automatically generate one based on the interior or exterior ZONE setting in the .rif file. In an interior zone, this default is a perspective view looking from the center of the bounding volume in the direction of the -X axis. To make better use of this first run, a plan-section view is generated instead. The view point is located above the model and the view target point directly beneath at a height selected to cut through as many apertures as possible to best illustrate the layout.

This is done so that the entire floor area, walls and most (usually all) of the windows are included in order that ambient values for those surfaces are sampled in the ambient file. Also, this image can also be used as a quick visual check of overall data quality and to spot any potential modelling anomalies.

#### Step 5: Creating Octree and Ambient Files

This step is important as it creates files required for the subsequent grid point analysis and \‘seeds\’ them with initial data relevant to that analysis. The ambient file is effectively a cache that is updated continuously as different views are rendered and more ambient information required (but limited to the original settings used to create it as described in the previous section). By generating a redundant interior view during this step, the ambient file will be generated with most of the information required for Step 6.

#### Step 6: Calculating Grid Point Values

In this step, the grid points are fed through rtrace to calculate their irradiance values. Whilst the rad command used in Step 5 requires only a small number of simplified settings, from which it automatically derives all the detailed parameter values it needs, the rtrace command requires a full set of detailed parameter values to be provided as arguments.

As mentioned previously, there appears to be only limited scope for customising parameter values between this and the previous step as the ambient file is sensitive only to the initial settings. However, by customising detailed parameter values early in Step 4 and including them in the .rif file together with an OPTFILE entry, it is possible to use exactly the same parameters in the rtrace command by simply including the .opt file in its command line. This is done by prefixing it with an @ sign, which causes Radiance to interpret the contents of the file as parameter values. This results in a command such as the following:

    rtrace -I -h -ov @PROJECT.opt PROJECT.oct < PROJECT.pts > PROJECT.dat


In this command, the -I option specifies the output as surface irradiance, the -h option specifies the output of header information and -ov specifies that the output should include just the irradiance value. Thus, it reads each grid point position and surface normal line-by-line from the PROJECT.pts file and outputs surface irradiance data for each point on a separate line in the PROJECT.dat file.

#### Step 7: Parsing Grid Point Values

The output from rtrace when using the -ov option actually generates the red, green and blue component surface irradiance values. To make these values directly comparable to the single reference horizontal surface irradiance give with the -B option in the gensky command, the three components are multiplied by the human eye\’s relative response to each colour band. The equation for this is as follows (Ward and Eydelberg-Vileshin, 2002):

$$Irradiance = (0.265 \times C_{red}) + (0.67 \times C_{green}) + (0.065 \times C_{blue}) \tag{3}$$

In the case of the particular model used here, internal surfaces are defined using non-colour-specific surface reflectances and windows using glazing transmittances, so all the internal surfaces are effectively shades of gray. The sky and ground glow RGB modifiers are all equal, making the sky and ground also shades of grey. Thus, each component value in the calculated results are the same so using this equation is mostly redundant as the modifiers all add up to 1.0 anyway. However, it is still used in order to accommodate future models with arbitrary surface properties and sky definitions.

The definition of daylight factor at a point is basically the ratio of total received illuminance/irradiance to the total available illuminance/irradiance from an unobstructed overcast sky, expressed as a percentage. Thus, by using the -B option in the gensky command and specifically setting it to a value of 100 to define the total unobstructed horizontal irradiance in W/m2, each grid point surface irradiance value is effectively already a daylight factor given as a percentage.

#### Step 8: Determining Correlation

A high potential correlation between the simplified GPU method and Radiance was immediately obvious early in the process when dynamically switching between the two sets of results within the prototype tool. The overall daylight distribution patterns matched very well, however there was an apparent linear offset between the two. Radiance results were consistently lower than those from the simplified GPU method, being around 75% of the absolute grid point values calculated by the simplified GPU method.

Thus, in addition to determining the correlation coefficient, the comparative analysis was also used to determine the nature of this linear offset, as well as if and why it varies with different parameters. This required deriving the following additional information from each run:

• The average daylight factor over the Radiance data grid,

• The average daylight factor over the GPU method data grid,

• The average and standard deviation of absolute differences between the values from the two methods at each grid point, and

• The root mean square error of the difference at each grid point.

As average absolute difference is not directly comparable between different grids, and even when daylight varies significantly across the same grid, an additional metric was calculated by dividing each absolute difference by the absolute daylight factor calculated by Radiance at each grid point, here called the relative difference. The resulting average relative difference allows for direct comparison between grids and across disparate values within the same grid.

# Results

The comparative results show that the correlation between raw daylight factor distributions generated by the simplified GPU method and Radiance using the same room configurations is surprisingly good, typically above 0.995 across a wide band of room parameter values.

The results of the two methods only begin to diverge under the following conditions:

• Correlation falls to zero when the work-plane is located above the head height of apertures. In the simplified GPU method there is no direct contribution from apertures located entirely below the work-plane, so the only value at each grid point is the average internally reflected component which makes the daylight distribution entirely uniform. In Radiance, diffuse light from apertures create areas of differential luminance/radiance across the ceiling, which affects the work-plane even when it is located above all apertures. Thus, even the slightest variation in the Radiance daylight distribution when compared to uniform values from the GPU method will result in a correlation of zero. However, these variations are trivially small and the average daylight factors across both data sets give similar results.

• Without using mkillum, correlation begins to fall below 0.995 when the total area of apertures becomes very small compared to the total internal surface area of the room. This is because of unevenness in Radiance results when ambient values alone are not sufficient to correctly model the internal daylight distribution. Under these conditions, even very high ambient settings and a calculation time of several hours does not produce sufficiently accurate results, as shown in Figure 5. However, replacing apertures with illum sources using pre-calculated luminance/radiance based on external sky conditions is one solution to this problem. To accommodate this, a switch from ambient calculations to the use of illum sources was implemented whenever the total aperture area is below 7.5% of the total internal surface area.

• Another condition which causes unevenness in Radiance results, and therefore a reduction in correlation, is when one dimension of the room is very small compared to the other two and the only windows are positioned in the smallest wall. For example, when the room is very long but also very thin, or when it is very large in plan but has a low ceiling height. This begins to occur when the smallest dimension is approximately 25% of the average of the other two.

The results of all parametric runs are summarised in Figure 15 and full details of each run with associated Radiance input and output files are available at: (URL removed to anonymise).

## Compensating for Relative Differences

Whilst the correlation is very good, the comparison of average relative differences between the daylight factor values from each method shows that raw Radiance results are in the order of 75% of the GPU method values. Closer investigation by isolating the direct and diffuse components in each method show that this difference occurs almost entirely in the diffuse component. When the internal reflectances for all materials are set to zero, the resulting direct-only raw sky components of the daylight factor in both methods match very closely. Figure 7 shows an example of just how closely direct-only raw value match across the entire grid.

This is seen more clearly in Figure 8 which shows the sectional profile of daylight factor values from both methods mapped against distance from an example window. The graph on the left shows raw daylight factors simulated on a 750mm work-plane and aligned with the window center, moving away from it in 250mm increments. The values both closest and furthest from the window generated by the GPU method are notably higher than Radiance values. However, setting the reflectance of all internal surfaces to zero and recalculating just the sky component, as shown in the graph on the right, shows a much tighter match.

A detailed inspection of the graphs in Figure 15 show that average relative difference values between the two methods vary primarily as a function of surface reflectance, window frame size and the percentage of window area beneath the work-plane. For all other parameters, when ambient issues in Radiance are discounted, the average relative difference is more consistent.

Looking at these three parameters across all aperture configurations, their change profiles in each case show the same characteristics. Thus, the first attempt to compensate for relative differences involved the generation of polynomials of best fit for each of the change profiles and then applying them as modifiers to the internally reflected component in the GPU method, both separately and together.

A brute force approach was then used to determine the effect of each modifier and their combination. This involved iterating through all the available run data across all aperture configurations and parameter values many times to find the modifier with the maximum overall correlation and minimum overall difference between grid points. Whilst a number of arrangements of these best fit functions did show a reduction in overall average absolute differences, they all resulted in an overall reduction in correlation.

After trying a range of different best-fit functions and multipliers, it was found that applying just a simple linear modifier of 0.52 to the internally reflected component of the GPU method achieved the highest overall correlation together with the minimum overall absolute difference between all grid point values across all aperture configurations and parameter values.

Obviously the true nature of differences in diffuse values between the two methods is very complex and, as the two methods calculate diffuse effects in very different ways, their absolute daylight factor values will never exactly match. However, by applying a simple linear modifier, the absolute values for most common room configurations with relatively small window frames, average surface reflectances and no single dimension less that 25% of the other two, can be made to match quite closely.

The effect of this linear modifier can be illustrated by revisiting the sectional profile graphs and overlaying the new modified GPU method values as a solid purple line, as shown in Figure 9. Now the two extreme cases closest and farthest from the window match much more closely and there are only mall differences in between. The sky component values remain unchanged as the modifier is only applied to the internally reflected component.

## Graphical Analysis

To gain insight into the results, a series of summarising graphs were generated for each parametric run and aperture configuration. The number of combinations is large, but a representative set of these graphs are included in Figure 14. Each summarising graph shows the following metrics:

• Correlation Coefficient (corrCoeff):
The correlation value between raw results from both methods as a blue line. This is an indicator of how well the overall daylight distribution patterns from the two methods match, with values closer to one (1.0) indicating a good match.

• Average Relative Ratio (diffRatio and modRatio):
This is the average daylight factor value given by the GPU method across all grid points divided by the average Radiance daylight factor value, to give the relative ratio of the two. Values closer to one (1.0) indicate a good match between absolute individual daylight factor values. The unmodified raw values are shown as a dotted red line (diffRatio) and values after applying the modifier as a solid red line (modRatio).

• Average Relative Difference (relDiff and modDiff):
This is the average absolute difference between the two daylight factor values across all grid points, normalised by dividing it by the average Radiance daylight factor value to give the relative difference. In this case, values closer to zero (0.0) indicate a good match between absolute individual daylight factor values. The unmodified raw values are shown as a dotted green line (relDiff) and values after applying the modifier as a solid red line (modDiff).

Before considering the entire result set, it is worth looking in more detail at some specific parameter runs in order to better illustrate what the summarised data means and how it was interpreted.

### The Effect of Window Size

The first example parameter to consider is shown in Figure 10 and details the effect of changing the window width from 250mm to 4750mm, in this case the full width of the wall in the +Y axis. The correlation over the whole work-plane grid between the GPU method and Radiance shows a steady line very near 1 for each window width. It also shows that the relative ratio of values is also steady at around 0.75 (red dotted line), except at window widths below 500mm. This is due to inadequate ambient sampling in Radiance resulting in splotchy and uneven results, even at very high settings.

The effect of applying the linear modifier of 0.52 to the internally reflected component is very clear in this example, showing that relative ratios move much closer to 1 and the relative difference closer to zero.

The effect of changing window height is very similar, though the impact of inadequate ambient sampling is more pronounced at heights below 500mm. Figure 11 summarises changes in height from 250mm up to the full height of the wall at 2250mm. In this example, the relative ratios and differences between the methods are slightly less steady, but still within 10% even at near full height.

### The Effect of Surface Reflectance

The second example parameter is important as it shows that, whilst correlation between the two methods remains very high, the relative ratio and difference values are significantly affected by high surface reflectances. This pattern of variation is consistent across all aperture configurations and highlights a limitation in how Equation 1 treats surface reflectance. As the room-averaged reflectance approaches 1, the calculated internally reflected component approaches infinity.

Much of this effect can actually be discounted as it is actually quite difficult to achieve surface reflectances of 0.85 and above with commonly used building materials. However, the graphs do show that there is a real effect between average reflectances of 0.5 and 0.85. Also of interest is how the modified and unmodified relative ratio and difference lines are affected by the linear modifier. This simply indicates that the overall effect of the internally reflected component is closely tied to surface reflectance, and the linear modifier reduces this effect.

### The Effect of Work Plane Height

The final example parameter to consider is the effect of work-plane height. As this increases above a window sill, each method treats the area of window below the work-plane very differently. In the GPU method, any area of window below the work-plane contributes no sky or externally reflected component and is only considered statistically as part of the internally reflected component. This continues as the work-plane is raised to be at or above the heads of all windows in the room, at which point there is no sky or externally reflected contribution to any part of the grid. As the only contribution is the room-averaged internally reflected component, the daylight grid will be entirely uniform.

Radiance, on the other hand, tracks diffuse radiation in and around the room to each surface, making it independent of work-plane height. Thus, even though the work-plane may be located above the heads of all windows in the room, there may be bright patches of luminance from the ceilings and high up the walls that influences daylight factor values.

Even though this may result in very small variations, the correlation will always be zero when compared to an entirely uniform grid. This is clearly shown in Figure 13 when the work-plane height reaches 2000mm, the correlation value falls to zero. Even at 1750mm, slightly below the window head height of 1800mm but able the frame, the correlation falls to 0.5.

Figure 14 shows what is actually happening here. The average daylight factors taken over both grids are very similar. In Radiance however, diffuse reflections from the ground and window sill illuminate an area of the ceiling above the window. As the work-plane is very high and close to the ceiling, it receives some illuminance from this bright patch which creates a slight variation across the Radiance grid. If we magnify the differences between each point on the two grids by 10, as shown on the right in Figure 14, the result is no correlation.

# Conclusion

The comparative analysis has shown that daylight factors across a horizontal work-plane, calculated using the GPU method implemented in this work and those from Radiance using the same simple room model, are very highly correlated. It has also shown how a small modifier applied to the internally reflected component of the GPU method can provide a significant reduction in overall differences between absolute daylight factor values across a wide parameter range when compared to Radiance simulations.

At the same time, this work has also verified some of the limitations of the split-flux method already noted by others, in that it is less suited to rooms with very high average surface reflectances and geometry that has one internal dimension less than 25% of the other two.

However, this method is significantly more stable than Radiance across a wider range of room, window and surface parameters, and when dealing with small or widely separated windows. Also, even though there are conditions under which the absolute values of daylight factors diverge, the high correlation between daylight distribution patterns generated by this method and Radiance suggests that the primary cause and effect relationships involved are being accurately captured and that changes to room configuration are being accurately reflected as relative changes in the distribution pattern.

As the method is also fast enough to allow for real-time visual feedback of daylight distribution as calculation parameters are dynamically manipulated by the user, the authors argue that this makes the prototype tool developed here a valuable educational resource - allowing users to interactively investigate the relationships between daylight distribution, room dimensions, surface properties, window sizes and their positions within the envelope.

This work has also shown that, when used with a relatively simple rectangular room model, the process of calculating spatial daylight distributions can be so highly optimised that it becomes a trivial component of animation frame updates. This provides significant potential for future development of the GPU method, and the prototype tool, to handle dynamic and cumulative sky luminance conditions, climate-based daylight modelling and more complex shading and glazing systems.

# References

2. Bourgeois, D., Reinhart, C.F., Ward G. (2008). A Standard Daylight Coefficient Model for Dynamic Daylighting Simulations. Building Research and Information 2008;36 (1):68e82.

3. BRE. (1986). Estimating Daylight in Buildings: Parts 1 and 2. BREPress.

4. Crisp, V.H.C., Littlefair, P.J. (1984). Average Daylight Factor Prediction, Proceedings of the CIBSE National Lighting Conference, Robinson College, Cambridge, UK.

5. Goral, C. M., Torrance, K. E., Greenberg, D. P., Battaile, B. (1984). Modeling the Interaction of Light Between Diffuse Surfaces, Computer Graphics, vol. 18, no. 3, pp. 213-222.

6. Hopkinson, R.G., Longmore, J., Petherbridge P. (1954). An Empirical Formula for the Computation of the Indirect Component of Daylight Factors. Trans. Illum. Eng. Soc. (London) 19, pp. 201.

7. Lynes, J.A. (1968). Principles of Natural Lighting. Applied Science Publishers, Ltd., London, pp. 129.

8. Lynes, J.A. (1979). A Sequence for Daylighting Design. Lighting Research & Technology, 11(2), pp. 102-106.

9. Mardeljevic, J. (1997). Validation of a Lighting Simulation Program: A Study Using Measured Sky Brightness Distributions, Proceedings of the 8th European Lighting Conference, Amsterdam, pp. 555-569.

10. Mardaljevic, J. (2000). Daylight simulation: validation, sky models and daylight coefficients. PhD Thesis. Leicester, UK: De Montfort University.

12. Osborne, J., Donn M. (2011). Defining parameters for a quality daylight simulation validation dataset. International Building Performance Simulation Association 2011. Sydney. Pp. 1481-1488.

13. Reinhart, C.F., LoVerso, V.R.M. (2010). A Rules of Thumb-based Design Sequence for Diffuse Daylight, Lighting Research and Technology 42(1), pp. 7-31.

14. TechPowerUp, (2017). GPU Database, Available online at https://www.techpowerup.com/gpudb/.

15. The U.S. DOE, (2017). EnergyPlus Engineering Reference: The Reference to EnergyPlus Calculations. Available online at https://energyplus.net/sites/default/files/pdfs_v8.3.0/EngineeringReference.pdf.

16. Versage, R., Melo, A.P., Lamberts, R. (2010). Impact of different daylighting simulation results on the prediction of total energy consumption. In Proceedings of the SimBuild 2010 Conference, New York, NY, USA, 11–13 August 2010; pp. 1–7.

17. Ward, G., Rubinstein, F. (1988). A New Technique for Computer Simulation of Illuminated Spaces, Journal of Illuminating Engineering Society, I, pp. 80-91.

18. Ward, G., Eydelberg-Vileshin E. (2002). Picture Perfect RGB Rendering Using Spectral Prefiltering and Sharp Color Primaries, Thirteenth Eurographics Workshop on Rendering; Debevec, P., Gibson, S. (Editors).

19. Whitted, T. (1979). An Improved Illumination Model for Shaded Display. Proceedings of the 6th annual conference on Computer graphics and interactive techniques.

20. Winkelmann, F., Selkowitz, S. (1985). Daylighting Simulation in the DOE-2 Building Energy Analysis Program. Energy and Buildings, v. l8, pp.271-286.

21. Yoon, Y.B., Lee, K.H. (2013). Comparing Radiosity and Split-flux illuminance algorithm in Energyplus. In Proceedings of the Biannual Autumn Conference of Korean Institute of Architectural Sustainable Environment and Building System (KIAEBS), Seoul, Korea.

22. Yoon, Y.B., Jeong, W.R., Lee, K.H. (2014). Window material daylighting performance assessment algorithm: Comparing radiosity and split-flux methods. Energies 2014, 7, pp. 2362–2376.

23. Yoon, Y.B., Manandhar R., Lee K.H. (2014). Comparative Study of Two Daylighting Analysis Methods with Regard to Window Orientation and Interior Wall Reflectance, Energies 2014, 7, pp. 5825-5846.