commit ee7fffb9a792e665176e510ee958697bff88c6a4
parent 842049974571f1a16d969396442abda46b2830bf
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 22 Jun 2018 10:21:14 +0200
Push further the volumic_power2 test
Diffstat:
2 files changed, 66 insertions(+), 33 deletions(-)
diff --git a/src/sdis_solve_Xd.h b/src/sdis_solve_Xd.h
@@ -218,8 +218,8 @@ XD(sample_reinjection_dir)
f2_normalize(dir, dir);
#else
/* Sample a random direction around the normal whose cosine is 1/sqrt(3). To
- * do so we sample a position onto a cone whose height is 1 and the radius of
- * its base is sqrt(2). */
+ * do so we sample a position onto a cone whose height is 1/sqrt(2) and the
+ * radius of its base is 1. */
float frame[9];
ASSERT(fX(is_normalized)(rwalk->hit.normal));
@@ -368,7 +368,7 @@ XD(trace_radiative_path)
res = RES_BAD_OP;
goto error;
}
- alpha = interface_side_get_specular_fraction(interf, &frag);
+ alpha = interface_side_get_specular_fraction(interf, &frag);
r = ssp_rng_canonical(rng);
if(r < alpha) { /* Sample specular part */
reflect_3d(dir, f3_minus(dir, dir), N);
@@ -500,13 +500,15 @@ XD(solid_solid_boundary_temperature)
delta_boundary_front = delta_front*sqrt(DIM);
delta_boundary_back = delta_back *sqrt(DIM);
- /* Sample a reinjection direction */
+ /* Sample a reinjection direction and reflect it around the normal. Then
+ * reflect them on the back side of the interfaces */
XD(sample_reinjection_dir)(rwalk, rng, dir0);
XD(reflect)(dir2, dir0, rwalk->hit.normal);
fX(minus)(dir1, dir0);
fX(minus)(dir3, dir2);
- /* Trace the dir0 and dir1 */
+ /* Trace the sampled directions to adjust the reinjection distance of the
+ * random walk on both sides of the interface. */
fX_set_dX(pos, rwalk->vtx.P);
f2(range0, 0, (float)delta_boundary_front*RAY_RANGE_MAX_SCALE);
f2(range1, 0, (float)delta_boundary_back *RAY_RANGE_MAX_SCALE);
diff --git a/src/test_sdis_volumic_power2.c b/src/test_sdis_volumic_power2.c
@@ -21,6 +21,7 @@
#define Pw 10000 /* Volumic power */
#define NONE -1
#define DELTA 0.01
+#define DELTA_PSQUARE 0.01
struct reference {
double pos[3];
@@ -28,6 +29,10 @@ struct reference {
double temperature_3d; /* In celcius */
};
+/* 100000 realisations
+ * (0 0.85); 2D = 190.29; 3D = 189.13 ~ 186.505 +/- 0.553322 [184.845, 188.165]
+ */
+
/* Results in Celcius with delta 0.01 and 10000 realisations
* 0.85: 2D = 190.29 ~ 194.022 +/- 1.87163 [188.407, 199.637]; #failures: 2
* 0.65: 2D = 259.95 ~ 265.535 +/- 2.14781 [259.091, 271.978]; #failures: 5
@@ -188,7 +193,6 @@ fluid_get_temperature
return ((const struct fluid*)sdis_data_cget(data))->temperature;
}
-
/*******************************************************************************
* Interfaces
******************************************************************************/
@@ -213,6 +217,40 @@ interface_get_temperature
return ((const struct interf*)sdis_data_cget(data))->temperature;
}
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+check(struct sdis_scene* scn, const struct reference refs[], const size_t nrefs)
+{
+ struct sdis_estimator* estimator = NULL;
+ struct sdis_mc T = SDIS_MC_NULL;
+ size_t nreals;
+ size_t nfails;
+ double pos[3] = {0,0};
+ size_t i;
+
+ FOR_EACH(i, 0, nrefs) {
+ double Tc;
+ pos[0] = refs[i].pos[0];
+ pos[1] = refs[i].pos[1];
+ pos[2] = refs[i].pos[2];
+
+ CHK(sdis_solve_probe(scn, N, pos, INF, 1.f, -1, 0, &estimator) == RES_OK);
+ CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
+ CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
+ CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
+ Tc = T.E - 273.15; /* Convert in Celcius */
+ printf("Temperature at (%g %g %g) = %g ~ %g +/- %g [%g, %g]\n",
+ SPLIT3(pos), refs[i].temperature_2d, Tc, T.SE, Tc-3*T.SE, Tc+3*T.SE);
+ printf("#realisations: %lu; #failures: %lu\n",
+ (unsigned long)nreals, (unsigned long)nfails);
+ /*CHK(eq_eps(Tc, refs[i].temperature, T.SE*3));*/
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+ }
+}
+
+
int
main(int argc, char** argv)
{
@@ -236,14 +274,8 @@ main(int argc, char** argv)
struct sdis_interface* interf_solid1_fluid1 = NULL;
struct sdis_interface* interf_solid1_fluid2 = NULL;
struct sdis_interface* interfaces[18 /*#rectangles*/];
- struct sdis_estimator* estimator = NULL;
- struct sdis_mc T = SDIS_MC_NULL;
- size_t nreals;
- size_t nfails;
- double pos[3] = {0,0,0};
- size_t i;
/* In celcius. Computed by EDF with Syrthes */
- const struct reference refs[] = { /* Lambda1=1, Lambda2=10, Pw = 10000 */
+ const struct reference refs1[] = { /* Lambda1=1, Lambda2=10, Pw = 10000 */
{{0, 0.85, 0}, 190.29, 189.13},
{{0, 0.65, 0}, 259.95, 247.09},
{{0, 0.45, 0}, 286.33, 308.42},
@@ -253,7 +285,12 @@ main(int argc, char** argv)
{{0,-0.35, 0}, 123.26, 123.43},
{{0,-0.55, 0}, 90.250, 90.040}
};
- size_t nrefs = sizeof(refs)/sizeof(struct reference);
+ const struct reference refs2[] = { /* Lambda1=0.1, Lambda2=10, Pw=10000 */
+ {{0, 0.85}, 678.170, -1},
+ {{0, 0.65}, 1520.84, -1},
+ {{0, 0.45}, 1794.57, -1},
+ {{0, 0.25}, 1429.74, -1}
+ };
(void)argc, (void)argv;
CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
@@ -309,7 +346,7 @@ main(int argc, char** argv)
solid_param->cp = 500000;
solid_param->rho = 1000;
solid_param->lambda = 10;
- solid_param->delta = DELTA;
+ solid_param->delta = DELTA_PSQUARE;
solid_param->P = Pw;
solid_param->T = -1;
CHK(sdis_solid_create(dev, &solid_shader, data, &solid2) == RES_OK);
@@ -388,25 +425,19 @@ main(int argc, char** argv)
exit(0);
#endif
- FOR_EACH(i, 0, nrefs) {
- double Tc;
- pos[0] = refs[i].pos[0];
- pos[1] = refs[i].pos[1];
- pos[2] = refs[i].pos[2];
+ printf(">>> Check 1\n");
+ check(scn, refs1, sizeof(refs1)/sizeof(struct reference));
- CHK(sdis_solve_probe(scn, N, pos, INF, 1.f, -1, 0, &estimator) == RES_OK);
- CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
- CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
- CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
- Tc = T.E - 273.15; /* Convert in Celcius */
- printf("Temperature at (%g %g); 2D = %g; 3D = %g ~ %g +/- %g [%g, %g]\n",
- SPLIT2(pos), refs[i].temperature_2d, refs[i].temperature_3d,
- Tc, T.SE, Tc-3*T.SE, Tc+3*T.SE);
- printf("#realisations: %lu; #failures: %lu\n",
- (unsigned long)nreals, (unsigned long)nfails);
- /*CHK(eq_eps(Tc, refs[i].temperature, T.SE*3));*/
- CHK(sdis_estimator_ref_put(estimator) == RES_OK);
- }
+ /* Update the scene */
+ CHK(sdis_scene_ref_put(scn) == RES_OK);
+ data = sdis_medium_get_data(solid1);
+ solid_param = sdis_data_get(data);
+ solid_param->lambda = 0.1;
+ CHK(sdis_scene_create(dev, ntriangles, get_indices, get_interface,
+ nvertices, get_position, interfaces, &scn) == RES_OK);
+
+ printf("\n>>> Check 2\n");
+ check(scn, refs2, sizeof(refs2)/sizeof(struct reference));
/* Release the interfaces */
CHK(sdis_interface_ref_put(interf_solid1_adiabatic) == RES_OK);
