test_ssol_solver6.c (10597B)
1 /* Copyright (C) 2018-2026 |Meso|Star> (contact@meso-star.com) 2 * Copyright (C) 2016, 2018 CNRS 3 * 4 * This program is free software: you can redistribute it and/or modify 5 * it under the terms of the GNU General Public License as published by 6 * the Free Software Foundation, either version 3 of the License, or 7 * (at your option) any later version. 8 * 9 * This program is distributed in the hope that it will be useful, 10 * but WITHOUT ANY WARRANTY; without even the implied warranty of 11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12 * GNU General Public License for more details. 13 * 14 * You should have received a copy of the GNU General Public License 15 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 16 17 #include "ssol.h" 18 #include "test_ssol_utils.h" 19 20 #define REFLECTIVITY 0 21 #include "test_ssol_materials.h" 22 23 #define PLANE_NAME SQUARE 24 #define HALF_X 5 25 #define HALF_Y 5 26 #include "test_ssol_rect_geometry.h" 27 28 #define POLYGON_NAME POLY 29 #define HALF_X 5 30 #define HALF_Y 5 31 #include "test_ssol_rect2D_geometry.h" 32 33 #include <rsys/double33.h> 34 35 #include <star/s3d.h> 36 #include <star/ssp.h> 37 38 static void 39 get_wlen(const size_t i, double* wlen, double* data, void* ctx) 40 { 41 double wavelengths[3] = { 1, 2, 3 }; 42 double intensities[3] = { 1, 0.8, 1 }; 43 CHK(i < 3); 44 (void)ctx; 45 *wlen = wavelengths[i]; 46 *data = intensities[i]; 47 } 48 49 int 50 main(int argc, char** argv) 51 { 52 struct mem_allocator allocator; 53 struct ssol_device* dev; 54 struct ssp_rng* rng; 55 const struct ssp_rng* rng_state; 56 struct ssol_object *m_object1; 57 struct ssol_object *m_object2; 58 struct ssol_object* t_object1; 59 struct ssol_object* t_object2; 60 struct ssol_instance* heliostat1; 61 struct ssol_instance* heliostat2; 62 struct ssol_instance* target1; 63 struct ssol_instance* target2; 64 struct ssol_scene* scene; 65 struct ssol_shape* square; 66 struct ssol_vertex_data attribs[1] = { SSOL_VERTEX_DATA_NULL__ }; 67 struct ssol_shape* quad_square; 68 struct ssol_carving carving = SSOL_CARVING_NULL; 69 struct ssol_quadric quadric = SSOL_QUADRIC_DEFAULT; 70 struct ssol_punched_surface punched = SSOL_PUNCHED_SURFACE_NULL; 71 struct ssol_material* m_mtl; 72 struct ssol_material* v_mtl; 73 struct ssol_material* bck_mtl; 74 struct ssol_mirror_shader m_shader = SSOL_MIRROR_SHADER_NULL; 75 struct ssol_matte_shader bck_shader = SSOL_MATTE_SHADER_NULL; 76 struct ssol_sun* sun; 77 struct ssol_spectrum* spectrum; 78 struct ssol_estimator* estimator; 79 struct ssol_estimator* estimator2; 80 struct ssol_mc_global mc_global; 81 struct ssol_mc_global mc_global2; 82 struct ssol_mc_receiver mc_rcv; 83 double dir[3]; 84 double transform[12]; /* 3x4 column major matrix */ 85 double sum_w, sum_w2, E, V, SE; 86 size_t count; 87 (void) argc, (void) argv; 88 89 mem_init_proxy_allocator(&allocator, &mem_default_allocator); 90 91 CHK(ssol_device_create 92 (NULL, &allocator, SSOL_NTHREADS_DEFAULT, 0, &dev) == RES_OK); 93 94 CHK(ssp_rng_create(&allocator, SSP_RNG_THREEFRY, &rng) == RES_OK); 95 CHK(ssol_spectrum_create(dev, &spectrum) == RES_OK); 96 CHK(ssol_spectrum_setup(spectrum, get_wlen, 3, NULL) == RES_OK); 97 CHK(ssol_sun_create_directional(dev, &sun) == RES_OK); 98 CHK(ssol_sun_set_direction(sun, d3(dir, 0, 0, -1)) == RES_OK); 99 CHK(ssol_sun_set_spectrum(sun, spectrum) == RES_OK); 100 CHK(ssol_sun_set_dni(sun, 1000) == RES_OK); 101 CHK(ssol_scene_create(dev, &scene) == RES_OK); 102 CHK(ssol_scene_attach_sun(scene, sun) == RES_OK); 103 104 /* Create scene content */ 105 106 CHK(ssol_shape_create_mesh(dev, &square) == RES_OK); 107 attribs[0].usage = SSOL_POSITION; 108 attribs[0].get = get_position; 109 CHK(ssol_mesh_setup(square, SQUARE_NTRIS__, get_ids, 110 SQUARE_NVERTS__, attribs, 1, (void*) &SQUARE_DESC__) == RES_OK); 111 112 CHK(ssol_shape_create_punched_surface(dev, &quad_square) == RES_OK); 113 carving.get = get_polygon_vertices; 114 carving.operation = SSOL_AND; 115 carving.nb_vertices = POLY_NVERTS__; 116 carving.context = &POLY_EDGES__; 117 quadric.type = SSOL_QUADRIC_PLANE; 118 punched.nb_carvings = 1; 119 punched.quadric = &quadric; 120 punched.carvings = &carving; 121 CHK(ssol_punched_surface_setup(quad_square, &punched) == RES_OK); 122 123 CHK(ssol_material_create_mirror(dev, &m_mtl) == RES_OK); 124 m_shader.normal = get_shader_normal; 125 m_shader.reflectivity = get_shader_reflectivity; 126 m_shader.roughness = get_shader_roughness; 127 CHK(ssol_mirror_setup(m_mtl, &m_shader, SSOL_MICROFACET_BECKMANN) == RES_OK); 128 CHK(ssol_material_create_matte(dev, &bck_mtl) == RES_OK); 129 bck_shader.normal = get_shader_normal; 130 bck_shader.reflectivity = get_shader_reflectivity_2; 131 CHK(ssol_matte_setup(bck_mtl, &bck_shader) == RES_OK); 132 CHK(ssol_material_create_virtual(dev, &v_mtl) == RES_OK); 133 134 /* 1st reflector */ 135 CHK(ssol_object_create(dev, &m_object1) == RES_OK); 136 CHK(ssol_object_add_shaded_shape(m_object1, quad_square, m_mtl, m_mtl) == RES_OK); 137 CHK(ssol_object_instantiate(m_object1, &heliostat1) == RES_OK); 138 CHK(ssol_scene_attach_instance(scene, heliostat1) == RES_OK); 139 d33_rotation_pitch(transform, PI); /* flip faces: invert normal */ 140 d3_splat(transform + 9, 0); 141 transform[9] = -25; 142 CHK(ssol_instance_set_transform(heliostat1, transform) == RES_OK); 143 144 /* 2nd reflector */ 145 CHK(ssol_object_create(dev, &m_object2) == RES_OK); 146 CHK(ssol_object_add_shaded_shape(m_object2, quad_square, m_mtl, m_mtl) == RES_OK); 147 CHK(ssol_object_instantiate(m_object2, &heliostat2) == RES_OK); 148 CHK(ssol_scene_attach_instance(scene, heliostat2) == RES_OK); 149 d33_rotation_pitch(transform, PI); /* flip faces: invert normal */ 150 d3_splat(transform + 9, 0); 151 transform[9] = +25; 152 CHK(ssol_instance_set_transform(heliostat2, transform) == RES_OK); 153 154 /* 1st target */ 155 CHK(ssol_object_create(dev, &t_object1) == RES_OK); 156 CHK(ssol_object_add_shaded_shape(t_object1, square, bck_mtl, v_mtl) == RES_OK); 157 CHK(ssol_object_instantiate(t_object1, &target1) == RES_OK); 158 CHK(ssol_instance_set_transform(target1, transform) == RES_OK); 159 CHK(ssol_instance_set_receiver(target1, SSOL_FRONT, 0) == RES_OK); 160 CHK(ssol_instance_sample(target1, 0) == RES_OK); 161 CHK(ssol_scene_attach_instance(scene, target1) == RES_OK); 162 d33_rotation_pitch(transform, PI); /* flip faces: invert normal */ 163 d3_splat(transform + 9, 0); 164 transform[9] = -25; 165 transform[11] = 10; 166 CHK(ssol_instance_set_transform(target1, transform) == RES_OK); 167 168 /* 2nd target */ 169 CHK(ssol_object_create(dev, &t_object2) == RES_OK); 170 CHK(ssol_object_add_shaded_shape(t_object2, square, bck_mtl, v_mtl) == RES_OK); 171 CHK(ssol_object_instantiate(t_object2, &target2) == RES_OK); 172 CHK(ssol_instance_set_transform(target2, transform) == RES_OK); 173 CHK(ssol_instance_set_receiver(target2, SSOL_BACK, 0) == RES_OK); 174 CHK(ssol_instance_sample(target2, 0) == RES_OK); 175 CHK(ssol_scene_attach_instance(scene, target2) == RES_OK); 176 d33_set_identity(transform); 177 d3_splat(transform + 9, 0); 178 transform[9] = +25; 179 transform[11] = 10; 180 CHK(ssol_instance_set_transform(target2, transform) == RES_OK); 181 182 #define N__ 10000 183 #define GET_MC_RCV ssol_estimator_get_mc_receiver 184 CHK(ssol_solve(scene, rng, N__, 0, NULL, &estimator) == RES_OK); 185 CHK(ssol_estimator_get_realisation_count(estimator, &count) == RES_OK); 186 CHK(count == N__); 187 CHK(ssol_estimator_get_mc_global(estimator, &mc_global) == RES_OK); 188 print_global(&mc_global); 189 CHK(eq_eps(mc_global.shadowed.E, 100000, 2 * 100000/sqrt(N__)) == 1); 190 CHK(eq_eps(mc_global.missing.E, 0, 0) == 1); 191 192 CHK(GET_MC_RCV(estimator, target1, SSOL_FRONT, &mc_rcv) == RES_OK); 193 printf("Ir(target1) = %g +/- %g\n", 194 mc_rcv.incoming_flux.E, mc_rcv.incoming_flux.SE); 195 CHK(eq_eps(mc_rcv.incoming_flux.E, 100000, 2*100000/sqrt(N__)) == 1); 196 CHK(mc_rcv.incoming_flux.E == mc_rcv.absorbed_flux.E); 197 198 CHK(GET_MC_RCV(estimator, target2, SSOL_BACK, &mc_rcv) == RES_OK); 199 printf("Ir(target2) = %g +/- %g\n", 200 mc_rcv.incoming_flux.E, mc_rcv.incoming_flux.SE); 201 CHK(eq_eps(mc_rcv.incoming_flux.E, 0, 1) == 1); 202 CHK(mc_rcv.incoming_flux.E == mc_rcv.absorbed_flux.E); 203 204 /* Launch another run that is statistically independent of the first one and 205 * that uses 3 more times samples */ 206 CHK(ssol_estimator_get_rng_state(estimator, &rng_state) == RES_OK); 207 CHK(ssol_solve(scene, rng_state, N__*3, 0, NULL, &estimator2) == RES_OK); 208 209 /* Check the estimator result */ 210 CHK(ssol_estimator_get_realisation_count(estimator2, &count) == RES_OK); 211 CHK(count == N__*3); 212 CHK(ssol_estimator_get_mc_global(estimator2, &mc_global2) == RES_OK); 213 CHK(eq_eps(mc_global2.shadowed.E, 100000, 2 * 100000/sqrt(N__)) == 1); 214 CHK(eq_eps(mc_global.shadowed.SE/sqrt(3), mc_global2.shadowed.SE, 1.e-1)); 215 216 CHK(mc_global.shadowed.E != mc_global2.shadowed.E); 217 CHK(mc_global.shadowed.SE != mc_global2.shadowed.SE); 218 219 /* Merge the 2 estimations */ 220 sum_w = mc_global.shadowed.E * N__; 221 sum_w += mc_global2.shadowed.E * N__*3; 222 V = mc_global.shadowed.SE * mc_global.shadowed.SE * N__; 223 sum_w2 = (V + mc_global.shadowed.E * mc_global.shadowed.E) * N__; 224 V = mc_global2.shadowed.SE * mc_global2.shadowed.SE * N__*3; 225 sum_w2 += (V + mc_global2.shadowed.E * mc_global2.shadowed.E) * N__*3; 226 E = sum_w / (4*N__); 227 V = sum_w2 / (4*N__) - E*E; 228 SE = sqrt(V/(4*N__)); 229 230 /* Check that the 2 runs are effectively independent, i.e. check the 231 * convergence ratio. By combining the 2 runs we have 4 more times samples, 232 * the standard deviation must be thus devided by 2 while the estimate must 233 * be compatible with the right value regarding the new error */ 234 CHK(eq_eps(E, 100000, 3*SE)); 235 CHK(eq_eps(SE, mc_global.shadowed.SE / 2, SE*1.e-2)); 236 237 /* Free data */ 238 CHK(ssol_instance_ref_put(heliostat1) == RES_OK); 239 CHK(ssol_instance_ref_put(target1) == RES_OK); 240 CHK(ssol_object_ref_put(m_object1) == RES_OK); 241 CHK(ssol_object_ref_put(t_object1) == RES_OK); 242 CHK(ssol_instance_ref_put(heliostat2) == RES_OK); 243 CHK(ssol_instance_ref_put(target2) == RES_OK); 244 CHK(ssol_object_ref_put(m_object2) == RES_OK); 245 CHK(ssol_object_ref_put(t_object2) == RES_OK); 246 CHK(ssol_shape_ref_put(square) == RES_OK); 247 CHK(ssol_shape_ref_put(quad_square) == RES_OK); 248 CHK(ssol_material_ref_put(m_mtl) == RES_OK); 249 CHK(ssol_material_ref_put(bck_mtl) == RES_OK); 250 CHK(ssol_material_ref_put(v_mtl) == RES_OK); 251 CHK(ssol_device_ref_put(dev) == RES_OK); 252 CHK(ssol_estimator_ref_put(estimator) == RES_OK); 253 CHK(ssol_estimator_ref_put(estimator2) == RES_OK); 254 CHK(ssol_scene_ref_put(scene) == RES_OK); 255 CHK(ssp_rng_ref_put(rng) == RES_OK); 256 CHK(ssol_spectrum_ref_put(spectrum) == RES_OK); 257 CHK(ssol_sun_ref_put(sun) == RES_OK); 258 259 check_memory_allocator(&allocator); 260 mem_shutdown_proxy_allocator(&allocator); 261 CHK(mem_allocated_size() == 0); 262 263 return 0; 264 }