commit d0cb8e268b33426274697dd4b119b5a824b22f41
parent e5617a4b72ef10ae98e59e5a8beee6a2bfe34c9d
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 29 May 2018 12:29:38 +0200
Implement a new 2D re-injection scheme of the random walk
Note that the 3D is disabled with since it is not correctly implemented
with this new scheme.
Diffstat:
3 files changed, 61 insertions(+), 51 deletions(-)
diff --git a/src/sdis_solve_Xd.h b/src/sdis_solve_Xd.h
@@ -377,7 +377,7 @@ XD(fluid_temperature)
return RES_OK;
}
-static void
+static INLINE void
XD(reinject_in_solid)
(const struct sdis_scene* scn,
const float fp_to_meter,
@@ -464,9 +464,7 @@ XD(solid_solid_boundary_temperature)
const struct sdis_medium* mdm;
struct s2d_hit hit0, hit1;
struct s2d_hit hit;
- double scale;
double power;
- double dst;
float range0[2], range1[2];
float dir0[2], dir1[2];
float pos[2];
@@ -474,24 +472,15 @@ XD(solid_solid_boundary_temperature)
/* Sample a direction */
r = ssp_rng_canonical(rng);
if(r < 0.5) {
- f3_set(dir0, rwalk->hit.normal);
- scale = 1.0;
- } else if(r < 0.75) {
dir0[0] = rwalk->hit.normal[0] - rwalk->hit.normal[1];
dir0[1] = rwalk->hit.normal[0] + rwalk->hit.normal[1];
- scale = sqrt(2.0);
} else {
dir0[0] = rwalk->hit.normal[0] + rwalk->hit.normal[1];
dir0[1] =-rwalk->hit.normal[0] + rwalk->hit.normal[1];
- scale = sqrt(2.0);
}
f2_normalize(dir0, dir0);
f2_minus(dir1, dir0);
- /* Adjust the delta */
- delta_front_boundary *= scale;
- delta_back_boundary *= scale;
-
/* Trace the dir0 and dir1 */
f2_set_d2(pos, rwalk->vtx.P);
f2(range0, 0, (float)delta_front_boundary*RAY_RANGE_MAX_SCALE);
@@ -500,13 +489,10 @@ XD(solid_solid_boundary_temperature)
S2D(scene_view_trace_ray(scn->s2d_view, pos, dir1, range1, &rwalk->hit, &hit1));
/* Define the reinjection distance */
- dst = MMIN
+ delta_boundary = MMIN
(MMIN(delta_front_boundary, delta_back_boundary),
MMIN(hit0.distance, hit1.distance));
- /* Define the delta to use */
- delta_boundary = dst / scale;
-
/* Define the reinjection side */
r = ssp_rng_canonical(rng);
proba = lambda_front / (lambda_front+lambda_back);
@@ -530,8 +516,8 @@ XD(solid_solid_boundary_temperature)
}
/* Reinject */
- XD(move_pos)(rwalk->vtx.P, dir, (float)dst);
- if(hit.distance == dst) {
+ XD(move_pos)(rwalk->vtx.P, dir, (float)delta_boundary);
+ if(hit.distance == delta_boundary) {
T->func = XD(boundary_temperature);
rwalk->mdm = NULL;
rwalk->hit = hit;
@@ -583,6 +569,7 @@ XD(solid_fluid_boundary_temperature)
const struct sdis_medium* mdm_back = NULL;
const struct sdis_medium* solid = NULL;
const struct sdis_medium* fluid = NULL;
+ struct sXd(hit) hit = SXD_HIT_NULL;
struct sdis_interface_fragment frag_fluid;
double hc;
double hr;
@@ -590,10 +577,13 @@ XD(solid_fluid_boundary_temperature)
double lambda;
double fluid_proba;
double radia_proba;
+ double delta;
double delta_boundary;
double r;
double tmp;
+ float pos[DIM];
float dir[DIM];
+ float range[2];
ASSERT(scn && fp_to_meter > 0 && rwalk && rng && T && ctx);
ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag));
@@ -617,7 +607,31 @@ XD(solid_fluid_boundary_temperature)
/* Fetch the solid properties */
lambda = solid_get_thermal_conductivity(solid, &rwalk->vtx);
- delta_boundary = solid_get_delta_boundary(solid, &rwalk->vtx);
+ delta = solid_get_delta(solid, &rwalk->vtx);
+ /* Note that elta boundary is *FIXED*. It MUST ensure that the orthogonal
+ * distance from the boundary to the point to chalenge is equal to delta. */
+ delta_boundary = sqrt(2.0) * delta;
+
+ /* Sample a direction. FIXME this only works in 2D */
+#if DIM == 3
+ FATAL("Un-implemented yet!\n");
+#endif
+ r = ssp_rng_canonical(rng);
+ if(r < 0.5) {
+ dir[0] = rwalk->hit.normal[0] - rwalk->hit.normal[1];
+ dir[1] = rwalk->hit.normal[0] + rwalk->hit.normal[1];
+ } else {
+ dir[0] = rwalk->hit.normal[0] + rwalk->hit.normal[1];
+ dir[1] =-rwalk->hit.normal[0] + rwalk->hit.normal[1];
+ }
+
+ /* Trace dir to adjust the reinection distance */
+ fX_set_dX(pos, rwalk->vtx.P);
+ f2(range, 0, (float)delta_boundary*RAY_RANGE_MAX_SCALE);
+ SXD(scene_view_trace_ray(scn->sXd(view), pos, dir, range, &rwalk->hit, &hit));
+
+ delta_boundary = MMIN(delta_boundary, hit.distance);
+ delta = delta_boundary / sqrt(2.0);
/* Fetch the boundary properties */
epsilon = interface_side_get_emissivity(interf, &frag_fluid);
@@ -627,7 +641,7 @@ XD(solid_fluid_boundary_temperature)
hr = 4.0 * BOLTZMANN_CONSTANT * ctx->Tref3 * epsilon;
/* Compute the probas to switch in solid or fluid random walk */
- tmp = lambda / (delta_boundary*fp_to_meter);
+ tmp = lambda / (delta*fp_to_meter);
fluid_proba = hc / (tmp + hr + hc);
radia_proba = hr / (tmp + hr + hc);
/*solid_proba = tmp / (tmp + hr + hc);*/
@@ -642,12 +656,27 @@ XD(solid_fluid_boundary_temperature)
rwalk->mdm = fluid;
rwalk->hit_side = rwalk->mdm == mdm_front ? SDIS_FRONT : SDIS_BACK;
} else { /* Solid random walk */
- rwalk->mdm = solid;
- fX(normalize)(dir, rwalk->hit.normal);
- if(solid == mdm_back) fX(minus)(dir, dir);
+ /* Handle the volumic power */
+ const double power = solid_get_volumic_power(solid, &rwalk->vtx);
+ if(power != SDIS_VOLUMIC_POWER_NONE) {
+ const double delta_in_meter = delta_boundary * fp_to_meter;
+ tmp = power * delta_in_meter * delta_in_meter / (2.0 * DIM * lambda);
+ T->value += tmp;
+ }
- /* "Reinject" the path into the solid along the surface normal. */
- XD(reinject_in_solid)(scn, 1, rwalk, delta_boundary, dir, T, 1);
+ /* Reinject */
+ XD(move_pos)(rwalk->vtx.P, dir, (float)delta_boundary);
+ if(hit.distance == delta_boundary) {
+ T->func = XD(boundary_temperature);
+ rwalk->mdm = NULL;
+ rwalk->hit = hit;
+ rwalk->hit_side = fX(dot)(hit.normal, dir) < 0 ? SDIS_FRONT : SDIS_BACK;
+ } else {
+ T->func = XD(solid_temperature);
+ rwalk->mdm = solid;
+ rwalk->hit = SXD_HIT_NULL;
+ rwalk->hit_side = SDIS_SIDE_NULL__;
+ }
}
}
@@ -744,20 +773,10 @@ XD(solid_temperature)
{
double position_start[DIM];
const struct sdis_medium* mdm;
- float low[DIM], upp[DIM];
- int i;
ASSERT(scn && fp_to_meter > 0 && rwalk && rng && T);
ASSERT(rwalk->mdm->type == SDIS_SOLID);
(void)ctx;
- /* FIXME hack */
- SXD(scene_view_get_aabb(scn->sXd(view), low, upp));
- FOR_EACH(i, 0, DIM) {
- low[i] *= low[i] < 0 ? 1.01f : 0.99f;
- upp[i] *= upp[i] < 0 ? 0.99f : 1.01f;
- }
-
-
/* Check the random walk consistency */
CHK(scene_get_medium(scn, rwalk->vtx.P, NULL, &mdm) == RES_OK);
if(mdm != rwalk->mdm) {
@@ -855,14 +874,6 @@ XD(solid_temperature)
/* Update the random walk position */
XD(move_pos)(rwalk->vtx.P, dir0, delta);
-#if 0
- FOR_EACH(i, 0, DIM) {
- if(rwalk->vtx.P[i] < low[i] || rwalk->vtx.P[i] > upp[i]) {
- log_err(scn->dev,"%s: invalid solid random walk.\n", FUNC_NAME);
- return RES_BAD_OP;
- }
- }
-#else
/* Fetch the current medium */
if(SXD_HIT_NONE(&rwalk->hit)) {
CHK(scene_get_medium(scn, rwalk->vtx.P, &info, &mdm) == RES_OK);
@@ -900,7 +911,6 @@ XD(solid_temperature)
#undef VEC_SPLIT
return RES_BAD_OP;
}
-#endif
/* Keep going while the solid random walk does not hit an interface */
} while(SXD_HIT_NONE(&rwalk->hit));
diff --git a/src/test_sdis_volumic_power2_2d.c b/src/test_sdis_volumic_power2_2d.c
@@ -17,9 +17,9 @@
#include "test_sdis_utils.h"
#include <rsys/math.h>
-#define N 10000 /* #realisations */
+#define N 40000 /* #realisations */
#define Pw 10000 /* Volumic power */
-#define MDb 1.0 /* Multiplier applied to delta to define delta_boundary */
+#define MDb sqrt(2.0) /* Multiplier applied to delta to define delta_boundary */
#define NONE -1
/* H delta T: expected 286.83 C */
diff --git a/src/test_sdis_volumic_power3_2d.c b/src/test_sdis_volumic_power3_2d.c
@@ -26,7 +26,7 @@
#define H1 5.0
#define H2 10.0
#define MDb 1.0
-#define N 10000 /* #realisations */
+#define N 400000 /* #realisations */
/*
* The 2D scene is composed of 3 stacked solid slabs whose middle slab has a
@@ -341,17 +341,17 @@ main(int argc, char** argv)
NULL, &data) == RES_OK);
interf_param = sdis_data_get(data);
interf_param->h = H2;
- interf_param->temperature = -1/*335.4141*/;
+ interf_param->temperature = 335.4141;
CHK(sdis_interface_create(dev, solid0, fluid, &interf_shader, data,
&interf_solid0_low) == RES_OK);
CHK(sdis_data_ref_put(data) == RES_OK);
- /* Create the solid0 T1 interace */
+ /* Create the solid0 upp interace */
CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
NULL, &data) == RES_OK);
interf_param = sdis_data_get(data);
interf_param->h = H1;
- interf_param->temperature = -1/*648.6217*/;
+ interf_param->temperature = 648.6217;
CHK(sdis_interface_create(dev, solid0, fluid, &interf_shader, data,
&interf_solid0_upp) == RES_OK);
CHK(sdis_data_ref_put(data) == RES_OK);
@@ -409,7 +409,7 @@ main(int argc, char** argv)
double Tref;
pos[0] = 0;
- pos[1] = 1.85 - (double)i*0.2;
+ pos[1] = 0.7; /*1.85 - (double)i*0.2;*/
ta = 1199.5651;
tb = 1207.1122;
