commit bbd08b5fb3c7b9998c30c38dce1de44efb312e70
parent ed2244e6f7032d0e252f4799c2781bd5598bdeee
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Thu, 21 Jan 2021 12:22:19 +0100
Improve radiative flux computation
Diffstat:
7 files changed, 98 insertions(+), 28 deletions(-)
diff --git a/src/sdis_Xd_begin.h b/src/sdis_Xd_begin.h
@@ -109,10 +109,11 @@ struct XD(temperature) {
struct XD(rwalk)* rwalk,
struct ssp_rng* rng,
struct XD(temperature)* temp);
+ void* ctx;
double value; /* Current value of the temperature */
int done;
};
-static const struct XD(temperature) XD(TEMPERATURE_NULL) = { NULL, 0, 0 };
+static const struct XD(temperature) XD(TEMPERATURE_NULL) = { NULL, NULL, 0, 0 };
#endif /* SDIX_<2|3>D_H */
diff --git a/src/sdis_heat_path_radiative_Xd.h b/src/sdis_heat_path_radiative_Xd.h
@@ -75,6 +75,7 @@ XD(trace_radiative_path)
if(SXD_HIT_NONE(&rwalk->hit)) { /* Fetch the ambient radiative temperature */
rwalk->hit_side = SDIS_SIDE_NULL__;
if(ctx->Tarad >= 0) {
+ struct radiative_path_ctx* rpctx = T->ctx;
T->value += ctx->Tarad;
T->done = 1;
@@ -94,6 +95,9 @@ XD(trace_radiative_path)
(ctx->heat_path, &vtx, T->value, SDIS_HEAT_VERTEX_RADIATIVE);
if(res != RES_OK) goto error;
}
+ if(rpctx && rpctx->status == FIRST_ABS_NOT_DEF_YET) {
+ rpctx->status = FIRST_ABS_OTHER;
+ }
break;
} else {
log_err(scn->dev,
@@ -140,8 +144,21 @@ XD(trace_radiative_path)
/* Switch in boundary temperature ? */
r = ssp_rng_canonical(rng);
if(r < epsilon) {
+ struct radiative_path_ctx* rpctx = T->ctx;
T->func = XD(boundary_path);
rwalk->mdm = NULL; /* The random walk is at an interface between 2 media */
+ if(rpctx && rpctx->status == FIRST_ABS_NOT_DEF_YET) {
+ /* Check if first absorption is with self or with other */
+ int is_self = rpctx->self
+ && htable_primitive_ids_find(rpctx->self, &rwalk->hit.prim.prim_id);
+ rpctx->status = is_self ? FIRST_ABS_SELF : FIRST_ABS_OTHER;
+ /* When computing radiative temperature and first absorption is with
+ * self, the end-of-path temperature is not used: just don't compute it! */
+ if(is_self) {
+ T->value = -1;
+ T->done = 1;
+ }
+ }
break;
}
diff --git a/src/sdis_misc.h b/src/sdis_misc.h
@@ -20,6 +20,35 @@
#include <rsys/float3.h>
#include <star/ssp.h>
+struct htable_primitive_ids;
+
+#include <rsys/hash_table.h>
+#define HTABLE_NAME primitive_ids
+#define HTABLE_KEY unsigned
+#define HTABLE_DATA char
+#include <rsys/hash_table.h>
+
+struct bound_flux_result {
+ double Tradiative;
+ double Tboundary;
+ double Tfluid;
+};
+#define BOUND_FLUX_RESULT_NULL__ {0,0,0}
+static const struct bound_flux_result
+BOUND_FLUX_RESULT_NULL = BOUND_FLUX_RESULT_NULL__;
+
+enum rad_path_first_abs {
+ FIRST_ABS_NOT_DEF_YET,
+ FIRST_ABS_SELF,
+ FIRST_ABS_OTHER
+};
+
+struct radiative_path_ctx {
+ struct htable_primitive_ids* self;
+ enum rad_path_first_abs status;
+};
+#define RADIATIVE_PATH_CTX_NULL__ { NULL, FIRST_ABS_NOT_DEF_YET }
+
struct accum {
double sum; /* Sum of MC weights */
double sum2; /* Sum of square MC weights */
diff --git a/src/sdis_realisation.h b/src/sdis_realisation.h
@@ -23,8 +23,11 @@
/* Forward declarations */
struct green_path_handle;
+struct sdis_heat_path;
struct sdis_scene;
struct ssp_rng;
+struct htable_primitive_ids;
+struct bound_flux_result;
enum flux_flag {
FLUX_FLAG_CONVECTIVE = BIT(FLUX_CONVECTIVE),
@@ -91,22 +94,24 @@ boundary_flux_realisation_2d
(struct sdis_scene* scn,
struct ssp_rng* rng,
const size_t iprim,
+ struct htable_primitive_ids* self,
const double uv[1],
const double time,
const enum sdis_side solid_side,
const int flux_mask, /* Combination of enum flux_flag */
- double weight[FLUX_NAMES_COUNT__]);
+ struct bound_flux_result* result);
extern LOCAL_SYM res_T
boundary_flux_realisation_3d
(struct sdis_scene* scn,
struct ssp_rng* rng,
const size_t iprim,
+ struct htable_primitive_ids* self,
const double uv[2],
const double time,
const enum sdis_side solid_side,
const int flux_mask, /* Combination of enum flux_flag */
- double weight[FLUX_NAMES_COUNT__]);
+ struct bound_flux_result* result);
/*******************************************************************************
* Realisation along a given ray at a given time. Available only in 3D.
diff --git a/src/sdis_realisation_Xd.h b/src/sdis_realisation_Xd.h
@@ -279,11 +279,12 @@ XD(boundary_flux_realisation)
(struct sdis_scene* scn,
struct ssp_rng* rng,
const size_t iprim,
+ struct htable_primitive_ids* self, /* NULL for probe computations */
const double uv[DIM],
const double time,
const enum sdis_side solid_side,
const int flux_mask,
- double weight[3])
+ struct bound_flux_result* result)
{
struct rwalk_context ctx = RWALK_CONTEXT_NULL;
struct XD(rwalk) rwalk;
@@ -292,6 +293,7 @@ XD(boundary_flux_realisation)
struct sXd(primitive) prim;
struct sdis_interface* interf = NULL;
struct sdis_medium* fluid_mdm = NULL;
+ struct radiative_path_ctx rpctx = RADIATIVE_PATH_CTX_NULL__;
#if SDIS_XD_DIMENSION == 2
float st;
@@ -307,7 +309,7 @@ XD(boundary_flux_realisation)
res_T res = RES_OK;
const char compute_radiative = (flux_mask & FLUX_FLAG_RADIATIVE) != 0;
const char compute_convective = (flux_mask & FLUX_FLAG_CONVECTIVE) != 0;
- ASSERT(uv && weight && time >= 0 );
+ ASSERT(uv && result && time >= 0 );
#if SDIS_XD_DIMENSION == 2
#define SET_PARAM(Dest, Src) (Dest).u = (Src);
@@ -341,6 +343,7 @@ XD(boundary_flux_realisation)
dX(set)(rwalk.vtx.P, P); \
fX(set)(rwalk.hit.normal, N); \
T = XD(TEMPERATURE_NULL); \
+ T.ctx = NULL; \
} (void)0
/* Compute boundary temperature */
@@ -348,7 +351,7 @@ XD(boundary_flux_realisation)
T.func = XD(boundary_path);
res = XD(compute_temperature)(scn, &ctx, &rwalk, rng, &T);
if(res != RES_OK) return res;
- weight[0] = T.value;
+ result->Tboundary = T.value;
/* Fetch the fluid medium */
interf = scene_get_interface(scn, (unsigned)iprim);
@@ -357,10 +360,18 @@ XD(boundary_flux_realisation)
/* Compute radiative temperature */
if(compute_radiative) {
RESET_WALK(fluid_side, fluid_mdm);
+ rpctx.self = self;
+ T.ctx = &rpctx;
T.func = XD(radiative_path);
res = XD(compute_temperature)(scn, &ctx, &rwalk, rng, &T);
if(res != RES_OK) return res;
- weight[1] = T.value;
+
+ if(rpctx.status == FIRST_ABS_SELF) {
+ result->Tradiative = -1;
+ } else {
+ ASSERT(T.value >= 0);
+ result->Tradiative = T.value;
+ }
}
/* Compute fluid temperature */
@@ -369,7 +380,7 @@ XD(boundary_flux_realisation)
T.func = XD(convective_path);
res = XD(compute_temperature)(scn, &ctx, &rwalk, rng, &T);
if(res != RES_OK) return res;
- weight[2] = T.value;
+ result->Tfluid = T.value;
}
#undef SET_PARAM
diff --git a/src/sdis_solve_boundary_Xd.h b/src/sdis_solve_boundary_Xd.h
@@ -446,6 +446,8 @@ XD(solve_boundary_flux)
int64_t irealisation;
size_t i;
size_t view_nprims;
+ struct htable_primitive_ids self;
+ int self_initialized = 0;
int progress = 0;
ATOMIC nsolved_realisations = 0;
ATOMIC res = RES_OK;
@@ -470,8 +472,12 @@ XD(solve_boundary_flux)
if(scene_is_2d(scn) != 0) { res = RES_BAD_ARG; goto error; }
#endif
+ htable_primitive_ids_init(scn->dev->allocator, &self);
+ self_initialized = 1;
SXD(scene_view_primitives_count(scn->sXd(view), &view_nprims));
FOR_EACH(i, 0, args->nprimitives) {
+ char one = 1;
+ unsigned prim;
if(args->primitives[i] >= view_nprims) {
log_err(scn->dev,
"%s: invalid primitive identifier `%lu'. It must be in the [0 %lu] range.\n",
@@ -481,6 +487,10 @@ XD(solve_boundary_flux)
res = RES_BAD_ARG;
goto error;
}
+ prim = (unsigned)args->primitives[i];
+ /* We don't reject multiple occurences */
+ res = htable_primitive_ids_set(&self, &prim, &one);
+ if(res != RES_OK) goto error;
}
/* Create the Star-XD shape of the boundary */
@@ -572,7 +582,7 @@ XD(solve_boundary_flux)
const struct sdis_interface* interf;
const struct sdis_medium *fmd, *bmd;
enum sdis_side solid_side, fluid_side;
- double T_brf[3] = { 0, 0, 0 };
+ struct bound_flux_result result = BOUND_FLUX_RESULT_NULL__;
const double Tref = scn->reference_temperature;
double epsilon, hc, hr, imposed_flux, imposed_temp;
size_t iprim;
@@ -655,8 +665,8 @@ XD(solve_boundary_flux)
flux_mask = 0;
if(hr > 0) flux_mask |= FLUX_FLAG_RADIATIVE;
if(hc > 0) flux_mask |= FLUX_FLAG_CONVECTIVE;
- res_simul = XD(boundary_flux_realisation)(scn, rng, iprim, uv, time,
- solid_side, flux_mask, T_brf);
+ res_simul = XD(boundary_flux_realisation)(scn, rng, iprim, &self, uv, time,
+ solid_side, flux_mask, &result);
/* Stop time registration */
time_sub(&t0, time_current(&t1), &t0);
@@ -666,16 +676,14 @@ XD(solve_boundary_flux)
continue;
} else if(res_simul == RES_OK) { /* Update accumulators */
const double usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
- const double Tboundary = T_brf[0];
- const double Tradiative = T_brf[1];
- const double Tfluid = T_brf[2];
- const double w_conv = hc * (Tboundary - Tfluid);
- const double w_rad = hr * (Tboundary - Tradiative);
+ const double w_conv = hc * (result.Tboundary - result.Tfluid);
+ const double w_rad = (result.Tradiative < 0) ?
+ 0 : hr * (result.Tboundary - result.Tradiative);
const double w_imp = (imposed_flux != SDIS_FLUX_NONE) ? imposed_flux : 0;
const double w_total = w_conv + w_rad + w_imp;
/* Temperature */
- acc_temp->sum += Tboundary;
- acc_temp->sum2 += Tboundary*Tboundary;
+ acc_temp->sum += result.Tboundary;
+ acc_temp->sum2 += result.Tboundary*result.Tboundary;
++acc_temp->count;
/* Time */
acc_time->sum += usec;
@@ -754,6 +762,7 @@ exit:
if(view) SXD(scene_view_ref_put(view));
if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
if(out_estimator) *out_estimator = estimator;
+ if(self_initialized) htable_primitive_ids_release(&self);
return (res_T)res;
error:
if(estimator) {
diff --git a/src/sdis_solve_probe_boundary_Xd.h b/src/sdis_solve_probe_boundary_Xd.h
@@ -479,7 +479,7 @@ XD(solve_probe_boundary_flux)
struct accum* acc_fimp = &acc_fi[ithread];
double time, epsilon, hc, hr, imposed_flux, imposed_temp;
int flux_mask = 0;
- double T_brf[3] = { 0, 0, 0 };
+ struct bound_flux_result result = BOUND_FLUX_RESULT_NULL__;
const double Tref = scn->reference_temperature;
size_t n;
int pcent;
@@ -513,8 +513,8 @@ XD(solve_probe_boundary_flux)
flux_mask = 0;
if(hr > 0) flux_mask |= FLUX_FLAG_RADIATIVE;
if(hc > 0) flux_mask |= FLUX_FLAG_CONVECTIVE;
- res_simul = XD(boundary_flux_realisation)(scn, rng, args->iprim, args->uv,
- time, solid_side, flux_mask, T_brf);
+ res_simul = XD(boundary_flux_realisation)(scn, rng, args->iprim, NULL,
+ args->uv, time, solid_side, flux_mask, &result);
/* Stop time registration */
time_sub(&t0, time_current(&t1), &t0);
@@ -524,16 +524,14 @@ XD(solve_probe_boundary_flux)
continue;
} else if(res_simul == RES_OK) { /* Update accumulators */
const double usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
- const double Tboundary = T_brf[0];
- const double Tradiative = T_brf[1];
- const double Tfluid = T_brf[2];
- const double w_conv = hc * (Tboundary - Tfluid);
- const double w_rad = hr * (Tboundary - Tradiative);
+ const double w_conv = hc * (result.Tboundary - result.Tfluid);
+ const double w_rad = (result.Tradiative < 0) ?
+ 0 : hr * (result.Tboundary - result.Tradiative);
const double w_imp = (imposed_flux != SDIS_FLUX_NONE) ? imposed_flux : 0;
const double w_total = w_conv + w_rad + w_imp;
/* Temperature */
- acc_temp->sum += Tboundary;
- acc_temp->sum2 += Tboundary*Tboundary;
+ acc_temp->sum += result.Tboundary;
+ acc_temp->sum2 += result.Tboundary*result.Tboundary;
++acc_temp->count;
/* Time */
acc_time->sum += usec;
