commit b82254f30264989868c8956574a0b5ad971fc0a6
parent 0db83e50e5bf68c1a4e301c354968bdb7e67332a
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Fri, 19 Oct 2018 13:47:06 +0200
Merge branch 'feature_flux_at_fluid_boundary' into develop
Diffstat:
10 files changed, 1227 insertions(+), 42 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -156,6 +156,7 @@ if(NOT NO_TEST)
new_test(test_sdis_solve_probe2_2d)
new_test(test_sdis_solve_probe3_2d)
new_test(test_sdis_solve_boundary)
+ new_test(test_sdis_solve_boundary_flux)
new_test(test_sdis_volumic_power)
# Additionnal tests
diff --git a/src/sdis.h b/src/sdis.h
@@ -76,6 +76,12 @@ enum sdis_medium_type {
SDIS_MEDIUM_TYPES_COUNT__
};
+enum sdis_estimator_type {
+ SDIS_TEMPERATURE_ESTIMATOR,
+ SDIS_FLUX_ESTIMATOR,
+ SDIS_EST_TYPES_COUNT__
+};
+
/* Random walk vertex, i.e. a spatiotemporal position at a given step of the
* random walk. */
struct sdis_rwalk_vertex {
@@ -336,7 +342,7 @@ SDIS_API res_T
sdis_accum_buffer_unmap
(const struct sdis_accum_buffer* buf);
-/* Helper function that matches the `sdis_write_accums_T' functor type. On can
+/* Helper function that matches the `sdis_write_accums_T' functor type. One can
* send this function directly to the sdis_solve_camera function, to fill the
* accum buffer with the estimation of the radiative temperature that reaches
* each pixel of an image whose definition matches the definition of the accum
@@ -527,9 +533,14 @@ sdis_estimator_ref_put
(struct sdis_estimator* estimator);
SDIS_API res_T
+sdis_estimator_get_type
+ (const struct sdis_estimator* estimator,
+ enum sdis_estimator_type* type);
+
+SDIS_API res_T
sdis_estimator_get_realisation_count
(const struct sdis_estimator* estimator,
- size_t* nrealisations);
+ size_t* nrealisations); /* Succesfull ones */
SDIS_API res_T
sdis_estimator_get_failure_count
@@ -541,6 +552,22 @@ sdis_estimator_get_temperature
(const struct sdis_estimator* estimator,
struct sdis_mc* temperature);
+SDIS_API res_T
+sdis_estimator_get_convective_flux
+ (const struct sdis_estimator* estimator,
+ struct sdis_mc* flux);
+
+SDIS_API res_T
+sdis_estimator_get_radiative_flux
+ (const struct sdis_estimator* estimator,
+ struct sdis_mc* flux);
+
+
+SDIS_API res_T
+sdis_estimator_get_total_flux
+ (const struct sdis_estimator* estimator,
+ struct sdis_mc* flux);
+
/*******************************************************************************
* Miscellaneous functions
******************************************************************************/
@@ -595,6 +622,31 @@ sdis_solve_boundary
const double reference_temperature, /* In Kelvin */
struct sdis_estimator** estimator);
+/* Flux solver */
+SDIS_API res_T
+sdis_solve_probe_boundary_flux
+ (struct sdis_scene* scn,
+ const size_t nrealisations, /* #realisations */
+ const size_t iprim, /* Identifier of the primitive on which the probe lies */
+ const double uv[2], /* Parametric coordinates of the probe onto the primitve */
+ const double time, /* Observation time */
+ const double fp_to_meter, /* Scale from floating point units to meters */
+ const double ambient_radiative_temperature, /* In Kelvin */
+ const double reference_temperature, /* In Kelvin */
+ struct sdis_estimator** estimator);
+
+SDIS_API res_T
+sdis_solve_boundary_flux
+ (struct sdis_scene* scn,
+ const size_t nrealisations, /* #realisations */
+ const size_t primitives[], /* List of boundary primitives to handle */
+ const size_t nprimitives, /* #primitives */
+ const double time, /* Observation time */
+ const double fp_to_meter, /* Scale from floating point units to meters */
+ const double ambient_radiative_temperature, /* In Kelvin */
+ const double reference_temperature, /* In Kelvin */
+ struct sdis_estimator** estimator);
+
END_DECLS
#endif /* SDIS_H */
diff --git a/src/sdis_estimator.c b/src/sdis_estimator.c
@@ -28,6 +28,8 @@ estimator_release(ref_T* ref)
ASSERT(ref);
estimator = CONTAINER_OF(ref, struct sdis_estimator, ref);
dev = estimator->dev;
+ ASSERT((estimator->fluxes!=NULL) == (estimator->type==SDIS_FLUX_ESTIMATOR));
+ MEM_RM(dev->allocator, estimator->fluxes);
MEM_RM(dev->allocator, estimator);
SDIS(device_ref_put(dev));
}
@@ -52,6 +54,15 @@ sdis_estimator_ref_put(struct sdis_estimator* estimator)
}
res_T
+sdis_estimator_get_type
+ (const struct sdis_estimator* estimator, enum sdis_estimator_type* type)
+{
+ if(!estimator || !type) return RES_BAD_ARG;
+ *type = estimator->type;
+ return RES_OK;
+}
+
+res_T
sdis_estimator_get_realisation_count
(const struct sdis_estimator* estimator, size_t* nrealisations)
{
@@ -78,16 +89,54 @@ sdis_estimator_get_temperature
return RES_OK;
}
+res_T
+sdis_estimator_get_convective_flux
+ (const struct sdis_estimator* estimator, struct sdis_mc* flux)
+{
+ if(!estimator || !flux ||estimator->type != SDIS_FLUX_ESTIMATOR)
+ return RES_BAD_ARG;
+ ASSERT(estimator->fluxes);
+ *flux = estimator->fluxes[FLUX_CONVECTIVE__];
+ return RES_OK;
+}
+
+res_T
+sdis_estimator_get_radiative_flux
+ (const struct sdis_estimator* estimator, struct sdis_mc* flux)
+{
+ if(!estimator || !flux || estimator->type != SDIS_FLUX_ESTIMATOR)
+ return RES_BAD_ARG;
+ ASSERT(estimator->fluxes);
+ *flux = estimator->fluxes[FLUX_RADIATIVE__];
+ return RES_OK;
+}
+
+res_T
+sdis_estimator_get_total_flux
+ (const struct sdis_estimator* estimator, struct sdis_mc* flux)
+{
+ if(!estimator || !flux || estimator->type != SDIS_FLUX_ESTIMATOR)
+ return RES_BAD_ARG;
+ ASSERT(estimator->fluxes);
+ *flux = estimator->fluxes[FLUX_TOTAL__];
+ return RES_OK;
+}
+
/*******************************************************************************
* Local functions
******************************************************************************/
res_T
-estimator_create(struct sdis_device* dev, struct sdis_estimator** out_estimator)
+estimator_create
+ (struct sdis_device* dev,
+ const enum sdis_estimator_type type,
+ struct sdis_estimator** out_estimator)
{
struct sdis_estimator* estimator = NULL;
res_T res = RES_OK;
- if(!dev || !out_estimator) {
+ if(!dev || !out_estimator
+ || (type != SDIS_TEMPERATURE_ESTIMATOR && type != SDIS_FLUX_ESTIMATOR))
+ {
res = RES_BAD_ARG;
goto error;
}
@@ -97,9 +146,13 @@ estimator_create(struct sdis_device* dev, struct sdis_estimator** out_estimator)
res = RES_MEM_ERR;
goto error;
}
+ estimator->type = type;
+ estimator->fluxes = (type != SDIS_FLUX_ESTIMATOR) ? NULL
+ : MEM_CALLOC(dev->allocator, FLUX_NAMES_COUNT__, sizeof(struct sdis_mc));
ref_init(&estimator->ref);
SDIS(device_ref_get(dev));
estimator->dev = dev;
+ if(type == SDIS_FLUX_ESTIMATOR && !estimator->fluxes) goto error;
exit:
if(out_estimator) *out_estimator = estimator;
diff --git a/src/sdis_estimator_c.h b/src/sdis_estimator_c.h
@@ -21,12 +21,22 @@
/* Forward declarations */
struct sdis_device;
struct sdis_estimator;
+enum sdis_estimator_type;
+
+enum flux_names {
+ FLUX_CONVECTIVE__,
+ FLUX_RADIATIVE__,
+ FLUX_TOTAL__,
+ FLUX_NAMES_COUNT__
+};
struct sdis_estimator {
struct sdis_mc temperature;
+ struct sdis_mc* fluxes;
size_t nrealisations;
size_t nfailures;
+ enum sdis_estimator_type type;
ref_T ref;
struct sdis_device* dev;
};
@@ -37,6 +47,7 @@ struct sdis_estimator {
extern LOCAL_SYM res_T
estimator_create
(struct sdis_device* dev,
+ const enum sdis_estimator_type type,
struct sdis_estimator** estimator);
#endif /* SDIS_PROBE_ESTIMATOR_C_H */
diff --git a/src/sdis_scene_Xd.h b/src/sdis_scene_Xd.h
@@ -585,11 +585,11 @@ XD(setup_enclosure_geometry)(struct sdis_scene* scn, struct sencXd(enclosure)* e
/* Compute the S/V ratio */
#if DIM == 2
- CALL(sXd(scene_view_compute_contour_length)(enc_data->sXd(view), &S));
- CALL(sXd(scene_view_compute_area)(enc_data->sXd(view), &V));
+ CALL(s2d_scene_view_compute_contour_length(enc_data->s2d_view, &S));
+ CALL(s2d_scene_view_compute_area(enc_data->s2d_view, &V));
#else
- CALL(sXd(scene_view_compute_area)(enc_data->sXd(view), &S));
- CALL(sXd(scene_view_compute_volume)(enc_data->sXd(view), &V));
+ CALL(s3d_scene_view_compute_area(enc_data->s3d_view, &S));
+ CALL(s3d_scene_view_compute_volume(enc_data->s3d_view, &V));
#endif
/* The volume of the enclosure is actually negative since Star-Enc ensures
* that the normal of its primitives point outward the enclosure. Take its
@@ -607,11 +607,11 @@ XD(setup_enclosure_geometry)(struct sdis_scene* scn, struct sencXd(enclosure)* e
if(res != RES_OK) goto error;
FOR_EACH(iprim, 0, nprims) {
#if DIM == 2
- SENCXD(enclosure_get_segment_global_id
- (enc, iprim, darray_uint_data_get(&enc_data->local2global)+iprim));
+ senc2d_enclosure_get_segment_global_id
+ (enc, iprim, darray_uint_data_get(&enc_data->local2global)+iprim);
#else
- SENCXD(enclosure_get_triangle_global_id
- (enc, iprim, darray_uint_data_get(&enc_data->local2global)+iprim));
+ senc_enclosure_get_triangle_global_id
+ (enc, iprim, darray_uint_data_get(&enc_data->local2global)+iprim);
#endif
}
diff --git a/src/sdis_solve.c b/src/sdis_solve.c
@@ -17,6 +17,7 @@
#include "sdis_camera.h"
#include "sdis_device_c.h"
#include "sdis_estimator_c.h"
+#include "sdis_interface_c.h"
#include "sdis_solve_Xd.h"
/* Generate the 2D solver */
@@ -64,7 +65,7 @@ solve_pixel
size_t N = 0; /* #realisations that do not fail */
size_t irealisation;
res_T res = RES_OK;
- ASSERT(scn && mdm && rng && cam && ipix && nrealisations);
+ ASSERT(scn && mdm && rng && cam && ipix && nrealisations && Tref >= 0);
ASSERT(pix_sz && pix_sz[0] > 0 && pix_sz[1] > 0);
FOR_EACH(irealisation, 0, nrealisations) {
@@ -125,7 +126,7 @@ solve_tile
size_t mcode; /* Morton code of the tile pixel */
size_t npixels;
res_T res = RES_OK;
- ASSERT(scn && rng && mdm && cam && spp && origin && accums);
+ ASSERT(scn && rng && mdm && cam && spp && origin && accums && Tref >= 0);
ASSERT(size &&size[0] && size[1]);
ASSERT(pix_sz && pix_sz[0] > 0 && pix_sz[1] > 0);
@@ -207,7 +208,7 @@ sdis_solve_probe
}
/* Create the estimator */
- res = estimator_create(scn->dev, &estimator);
+ res = estimator_create(scn->dev, SDIS_TEMPERATURE_ESTIMATOR, &estimator);
if(res != RES_OK) goto error;
/* Retrieve the medium in which the submitted position lies */
@@ -299,10 +300,10 @@ sdis_solve_probe_boundary
/* Check the primitive identifier */
if(iprim >= scene_get_primitives_count(scn)) {
log_err(scn->dev,
-"%s: invalid primitive identifier `%lu'. It must be less than %lu.\n",
+"%s: invalid primitive identifier `%lu'. It must be in the [0 %lu] range.\n",
FUNC_NAME,
(unsigned long)iprim,
- (unsigned long)scene_get_primitives_count(scn));
+ (unsigned long)scene_get_primitives_count(scn)-1);
res = RES_BAD_ARG;
goto error;
}
@@ -349,7 +350,7 @@ sdis_solve_probe_boundary
}
/* Create the estimator */
- res = estimator_create(scn->dev, &estimator);
+ res = estimator_create(scn->dev, SDIS_TEMPERATURE_ESTIMATOR, &estimator);
if(res != RES_OK) goto error;
/* Here we go! Launch the Monte Carlo estimation */
@@ -566,3 +567,209 @@ sdis_solve_boundary
return res;
}
+res_T
+sdis_solve_probe_boundary_flux
+ (struct sdis_scene* scn,
+ const size_t nrealisations, /* #realisations */
+ const size_t iprim, /* Identifier of the primitive on which the probe lies */
+ const double uv[2], /* Parametric coordinates of the probe onto the primitve */
+ const double time, /* Observation time */
+ const double fp_to_meter, /* Scale from floating point units to meters */
+ const double Tarad, /* In Kelvin */
+ const double Tref, /* In Kelvin */
+ struct sdis_estimator** out_estimator)
+{
+ struct sdis_estimator* estimator = NULL;
+ struct ssp_rng_proxy* rng_proxy = NULL;
+ struct ssp_rng** rngs = NULL;
+ const struct sdis_interface* interf;
+ const struct sdis_medium *fmd, *bmd;
+ enum sdis_side solid_side, fluid_side;
+ double weight_t = 0, sqr_weight_t = 0;
+ double weight_fc = 0, sqr_weight_fc = 0;
+ double weight_fr = 0, sqr_weight_fr = 0;
+ double weight_f= 0, sqr_weight_f = 0;
+ double epsilon, hc, hr;
+ const int64_t rcount = (int64_t)nrealisations;
+ int64_t irealisation = 0;
+ size_t N = 0; /* #realisations that do not fail */
+ size_t i;
+ ATOMIC res = RES_OK;
+
+ if(!scn || !nrealisations || nrealisations > INT64_MAX || !uv || time < 0
+ || fp_to_meter <= 0 || Tref < 0
+ || !out_estimator) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Check the primitive identifier */
+ if(iprim >= scene_get_primitives_count(scn)) {
+ log_err(scn->dev,
+ "%s: invalid primitive identifier `%lu'. It must be in the [0 %lu] range.\n",
+ FUNC_NAME,
+ (unsigned long)iprim,
+ (unsigned long)scene_get_primitives_count(scn)-1);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Check parametric coordinates */
+ if(scene_is_2d(scn)) {
+ const double v = CLAMP(1.0 - uv[0], 0, 1);
+ if(uv[0] < 0 || uv[0] > 1 || !eq_eps(uv[0] + v, 1, 1.e-6)) {
+ log_err(scn->dev,
+ "%s: invalid parametric coordinates %g.\n"
+ "u + (1-u) must be equal to 1 with u [0, 1].\n",
+ FUNC_NAME, uv[0]);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ } else {
+ const double w = CLAMP(1 - uv[0] - uv[1], 0, 1);
+ if(uv[0] < 0 || uv[1] < 0 || uv[0] > 1 || uv[1] > 1
+ || !eq_eps(w + uv[0] + uv[1], 1, 1.e-6)) {
+ log_err(scn->dev,
+ "%s: invalid parametric coordinates [%g, %g].\n"
+ "u + v + (1-u-v) must be equal to 1 with u and v in [0, 1].\n",
+ FUNC_NAME, uv[0], uv[1]);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ }
+ /* Check medium is fluid on one side and solid on the other */
+ interf = scene_get_interface(scn, (unsigned long)iprim);
+ fmd = interface_get_medium(interf, SDIS_FRONT);
+ bmd = interface_get_medium(interf, SDIS_BACK);
+ if(!fmd || !bmd
+ || (!(fmd->type == SDIS_FLUID && bmd->type == SDIS_SOLID)
+ && !(fmd->type == SDIS_SOLID && bmd->type == SDIS_FLUID)))
+ {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ solid_side = (fmd->type == SDIS_SOLID) ? SDIS_FRONT : SDIS_BACK;
+ fluid_side = (fmd->type == SDIS_FLUID) ? SDIS_FRONT : SDIS_BACK;
+
+ /* Create the proxy RNG */
+ res = ssp_rng_proxy_create(scn->dev->allocator, &ssp_rng_mt19937_64,
+ scn->dev->nthreads, &rng_proxy);
+ if(res != RES_OK) goto error;
+
+ /* Create the per thread RNG */
+ rngs = MEM_CALLOC
+ (scn->dev->allocator, scn->dev->nthreads, sizeof(struct ssp_rng*));
+ if(!rngs) {
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ FOR_EACH(i, 0, scn->dev->nthreads) {
+ res = ssp_rng_proxy_create_rng(rng_proxy, i, rngs + i);
+ if(res != RES_OK) goto error;
+ }
+
+ /* Compute hr and hc */
+ if(scene_is_2d(scn)) {
+ res = interface_get_hc_epsilon_2d(&hc, &epsilon, scn, (unsigned long)iprim,
+ uv, time, fluid_side);
+ } else {
+ res = interface_get_hc_epsilon_3d(&hc, &epsilon, scn, (unsigned long)iprim,
+ uv, time, fluid_side);
+ }
+ hr = 4.0 * BOLTZMANN_CONSTANT * Tref * Tref * Tref * epsilon;
+
+ /* Create the estimator */
+ res = estimator_create(scn->dev, SDIS_FLUX_ESTIMATOR, &estimator);
+ if(res != RES_OK) goto error;
+
+ /* Here we go! Launch the Monte Carlo estimation */
+ omp_set_num_threads((int)scn->dev->nthreads);
+ #pragma omp parallel for schedule(static) reduction(+:weight_t,sqr_weight_t,\
+ weight_fc,sqr_weight_fc,weight_fr,sqr_weight_fr,weight_f,sqr_weight_f,N)
+ for(irealisation = 0; irealisation < rcount; ++irealisation) {
+ res_T res_local;
+ double T_brf[3] = { 0, 0, 0 };
+ const int ithread = omp_get_thread_num();
+ struct ssp_rng* rng = rngs[ithread];
+
+ if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occurred */
+
+ /* Fluid, Radiative and Solid temperatures */
+ if(scene_is_2d(scn)) {
+ res_local = probe_flux_realisation_2d(scn, rng, iprim, uv, time,
+ solid_side, fp_to_meter, Tarad, Tref, hr>0, hc>0, T_brf);
+ } else {
+ res_local = probe_flux_realisation_3d(scn, rng, iprim, uv, time,
+ solid_side, fp_to_meter, Tarad, Tref, hr>0, hc>0, T_brf);
+ }
+ if(res_local != RES_OK) {
+ if(res_local != RES_BAD_OP) {
+ ATOMIC_SET(&res, res_local);
+ continue;
+ }
+ } else {
+ 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_total = w_conv + w_rad;
+ weight_t += Tboundary;
+ sqr_weight_t += Tboundary * Tboundary;
+ weight_fc += w_conv;
+ sqr_weight_fc += w_conv * w_conv;
+ weight_fr += w_rad;
+ sqr_weight_fr += w_rad * w_rad;
+ weight_f += w_total;
+ sqr_weight_f += w_total * w_total;
+ ++N;
+ }
+ }
+ if(res != RES_OK) goto error;
+
+ setup_estimator(estimator, nrealisations, N, weight_t, sqr_weight_t);
+ setup_estimator_flux(estimator, FLUX_CONVECTIVE__, weight_fc, sqr_weight_fc);
+ setup_estimator_flux(estimator, FLUX_RADIATIVE__, weight_fr, sqr_weight_fr);
+ setup_estimator_flux(estimator, FLUX_TOTAL__, weight_f, sqr_weight_f);
+
+exit:
+ if(rngs) {
+ FOR_EACH(i, 0, scn->dev->nthreads) {
+ if(rngs[i]) SSP(rng_ref_put(rngs[i]));
+ }
+ MEM_RM(scn->dev->allocator, rngs);
+ }
+ if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
+ if(out_estimator) *out_estimator = estimator;
+ return (res_T)res;
+error:
+ if(estimator) {
+ SDIS(estimator_ref_put(estimator));
+ estimator = NULL;
+ }
+ goto exit;
+}
+
+res_T
+sdis_solve_boundary_flux
+ (struct sdis_scene* scn,
+ const size_t nrealisations, /* #realisations */
+ const size_t primitives [], /* List of boundary primitives to handle */
+ const size_t nprimitives, /* #primitives */
+ const double time, /* Observation time */
+ const double fp_to_meter, /* Scale from floating point units to meters */
+ const double Tarad, /* In Kelvin */
+ const double Tref, /* In Kelvin */
+ struct sdis_estimator** out_estimator)
+{
+ res_T res = RES_OK;
+ if(!scn) return RES_BAD_ARG;
+ if(scene_is_2d(scn)) {
+ res = solve_boundary_flux_2d(scn, nrealisations, primitives, nprimitives,
+ time, fp_to_meter, Tarad, Tref, out_estimator);
+ } else {
+ res = solve_boundary_flux_3d(scn, nrealisations, primitives, nprimitives,
+ time, fp_to_meter, Tarad, Tref, out_estimator);
+ }
+ return res;
+}
diff --git a/src/sdis_solve_Xd.h b/src/sdis_solve_Xd.h
@@ -97,7 +97,26 @@ setup_estimator
accum_sqr_weights / (double)nsuccesses
- estimator->temperature.E * estimator->temperature.E;
estimator->temperature.V = MMAX(estimator->temperature.V, 0);
- estimator->temperature.SE = sqrt(estimator->temperature.V / (double)nsuccesses);
+ estimator->temperature.SE =
+ sqrt(estimator->temperature.V / (double)nsuccesses);
+}
+
+static INLINE void
+setup_estimator_flux
+ (struct sdis_estimator* estimator,
+ const enum flux_names name,
+ const double accum_weights,
+ const double accum_sqr_weights)
+{
+ ASSERT(estimator && (unsigned)name < FLUX_NAMES_COUNT__ && estimator->fluxes);
+ ASSERT(estimator->nrealisations);
+ estimator->fluxes[name].E = accum_weights / (double)estimator->nrealisations;
+ estimator->fluxes[name].V =
+ accum_sqr_weights / (double)estimator->nrealisations
+ - estimator->fluxes[name].E * estimator->fluxes[name].E;
+ estimator->fluxes[name].V = MMAX(estimator->fluxes[name].V, 0);
+ estimator->fluxes[name].SE =
+ sqrt(estimator->fluxes[name].V / (double)estimator->nrealisations);
}
#endif /* SDIS_SOLVE_XD_H */
@@ -126,6 +145,7 @@ setup_estimator
#define SXD_HIT_NULL__ CONCAT(CONCAT(S, DIM), D_HIT_NULL__)
#define SXD_POSITION CONCAT(CONCAT(S, DIM), D_POSITION)
#define SXD_GEOMETRY_NORMAL CONCAT(CONCAT(S, DIM), D_GEOMETRY_NORMAL)
+#define SXD_GEOMETRY_NORMAL CONCAT(CONCAT(S, DIM), D_GEOMETRY_NORMAL)
#define SXD_VERTEX_DATA_NULL CONCAT(CONCAT(S, DIM), D_VERTEX_DATA_NULL)
#define SXD CONCAT(CONCAT(S, DIM), D)
#define SXD_FLOAT2 CONCAT(CONCAT(S, DIM), D_FLOAT2)
@@ -236,11 +256,11 @@ XD(boundary_get_position)(const unsigned ivert, float pos[DIM], void* context)
iprim = (unsigned)ctx->primitives[iprim_id];
SXD(scene_view_get_primitive(ctx->view, iprim, &prim));
#if DIM == 2
- SXD(segment_get_vertex_attrib(&prim, iprim_vert, SXD_POSITION, &attr));
- ASSERT(attr.type == SXD_FLOAT2);
+ s2d_segment_get_vertex_attrib(&prim, iprim_vert, S2D_POSITION, &attr);
+ ASSERT(attr.type == S2D_FLOAT2);
#else
- SXD(triangle_get_vertex_attrib(&prim, iprim_vert, SXD_POSITION, &attr));
- ASSERT(attr.type == SXD_FLOAT3);
+ s3d_triangle_get_vertex_attrib(&prim, iprim_vert, S3D_POSITION, &attr);
+ ASSERT(attr.type == S3D_FLOAT3);
#endif
fX(set)(pos, attr.value);
}
@@ -1273,7 +1293,7 @@ XD(solid_temperature)
tmp += (power*delta_s_in_meter*delta_s_in_meter)/(6*lambda) * tmp1;
#endif
- } else if (h == delta_solid) {
+ } else if(h == delta_solid) {
tmp += -(delta_s_in_meter*delta_s_in_meter*power)/(2.0*DIM*lambda);
}
T->value += tmp;
@@ -1461,7 +1481,7 @@ XD(boundary_realisation)
(struct sdis_scene* scn,
struct ssp_rng* rng,
const size_t iprim,
- const double uv[DIM],
+ const double uv[2],
const double time,
const enum sdis_side side,
const double fp_to_meter,
@@ -1479,7 +1499,7 @@ XD(boundary_realisation)
float st[2];
#endif
res_T res = RES_OK;
- ASSERT(uv && fp_to_meter > 0 && weight && time >= 0);
+ ASSERT(uv && fp_to_meter > 0 && weight && time >= 0 && Tref >= 0);
T.func = XD(boundary_temperature);
@@ -1522,6 +1542,163 @@ XD(boundary_realisation)
return RES_OK;
}
+static res_T
+XD(probe_flux_realisation)
+ (struct sdis_scene* scn,
+ struct ssp_rng* rng,
+ const size_t iprim,
+ const double uv[DIM],
+ const double time,
+ const enum sdis_side solid_side,
+ const double fp_to_meter,
+ const double Tarad,
+ const double Tref,
+ const char compute_radiative,
+ const char compute_convective,
+ double weight[3])
+{
+ struct rwalk_context ctx;
+ struct XD(rwalk) rwalk;
+ struct XD(temperature) T;
+ struct sXd(attrib) attr;
+ struct sXd(primitive) prim;
+#if SDIS_SOLVE_DIMENSION == 2
+ float st;
+#else
+ float st[2];
+#endif
+ double P[SDIS_SOLVE_DIMENSION];
+ float N[SDIS_SOLVE_DIMENSION];
+ const double Tr3 = Tref * Tref * Tref;
+ const enum sdis_side fluid_side =
+ (solid_side == SDIS_FRONT) ? SDIS_BACK : SDIS_FRONT;
+ res_T res = RES_OK;
+ ASSERT(uv && fp_to_meter > 0 && weight && time >= 0 && Tref >= 0);
+
+#if SDIS_SOLVE_DIMENSION == 2
+ #define SET_PARAM(Dest, Src) (Dest).u = (Src);
+ st = (float)uv[0];
+#else
+ #define SET_PARAM(Dest, Src) f2_set((Dest).uv, (Src));
+ f2_set_d2(st, uv);
+#endif
+
+ /* Fetch the primitive */
+ SXD(scene_view_get_primitive(scn->sXd(view), (unsigned int)iprim, &prim));
+
+ /* Retrieve the world space position of the probe onto the primitive */
+ SXD(primitive_get_attrib(&prim, SXD_POSITION, st, &attr));
+ dX_set_fX(P, attr.value);
+
+ /* Retrieve the primitive normal */
+ SXD(primitive_get_attrib(&prim, SXD_GEOMETRY_NORMAL, st, &attr));
+ fX(set)(N, attr.value);
+
+ #define RESET_WALK(Side, Mdm) \
+ rwalk = XD(RWALK_NULL); \
+ rwalk.hit_side = (Side); \
+ rwalk.hit.distance = 0; \
+ rwalk.vtx.time = time; \
+ rwalk.mdm = (Mdm); \
+ SET_PARAM(rwalk.hit, st); \
+ ctx.Tarad = Tarad; \
+ ctx.Tref3 = Tr3; \
+ rwalk.hit.prim = prim; \
+ dX(set)(rwalk.vtx.P, P); \
+ fX(set)(rwalk.hit.normal, N); \
+ T = XD(TEMPERATURE_NULL);
+
+ /* Compute boundary temperature */
+ RESET_WALK(solid_side, NULL);
+ T.func = XD(boundary_temperature);
+ res = XD(compute_temperature)(scn, fp_to_meter, &ctx, &rwalk, rng, &T);
+ if(res != RES_OK) return res;
+ weight[0] = T.value;
+
+ /* Compute radiative temperature */
+ if(compute_radiative) {
+ RESET_WALK(fluid_side, NULL);
+ T.func = XD(radiative_temperature);
+ res = XD(compute_temperature)(scn, fp_to_meter, &ctx, &rwalk, rng, &T);
+ if(res != RES_OK) return res;
+ weight[1] = T.value;
+ }
+
+ /* Compute fluid temperature */
+ if(compute_convective) {
+ const struct sdis_interface* interf =
+ scene_get_interface(scn, (unsigned long)iprim);
+ const struct sdis_medium* mdm = interface_get_medium(interf, fluid_side);
+
+ RESET_WALK(fluid_side, mdm);
+ T.func = XD(fluid_temperature);
+ res = XD(compute_temperature)(scn, fp_to_meter, &ctx, &rwalk, rng, &T);
+ if(res != RES_OK) return res;
+ weight[2] = T.value;
+ }
+
+ #undef SET_PARAM
+ #undef RESET_WALK
+
+ return RES_OK;
+}
+
+static res_T
+XD(interface_get_hc_epsilon)
+ (double *hc,
+ double* epsilon,
+ const struct sdis_scene* scn,
+ const unsigned iprim,
+ const double* uv,
+ const double time,
+ const enum sdis_side fluid_side)
+{
+ struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
+ struct sXd(attrib) attr;
+ struct sXd(primitive) prim;
+ struct sXd(hit) hit;
+ struct sdis_rwalk_vertex vtx;
+ const struct sdis_interface* interf;
+#if SDIS_SOLVE_DIMENSION == 2
+ float st;
+#else
+ float st[2];
+#endif
+ res_T res = RES_OK;
+
+ ASSERT(fluid_side == SDIS_FRONT || fluid_side == SDIS_BACK);
+
+#if SDIS_SOLVE_DIMENSION == 2
+ #define SET_PARAM(Dest, Src) (Dest).u = (Src);
+ st = (float) uv[0];
+#else
+ #define SET_PARAM(Dest, Src) f2_set((Dest).uv, (Src));
+ f2_set_d2(st, uv);
+#endif
+ res = sXd(scene_view_get_primitive(scn->sXd(view), iprim, &prim));
+ if(res != RES_OK) return res;
+ res = sXd(primitive_get_attrib(&prim, SXD_POSITION, st, &attr));
+ if(res != RES_OK) return res;
+ dX_set_fX(vtx.P, attr.value);
+ res = sXd(primitive_get_attrib(&prim, SXD_GEOMETRY_NORMAL, st, &attr));
+ if(res != RES_OK) return res;
+ fX(set)(hit.normal, attr.value);
+
+ hit.distance = 0;
+ hit.normal;
+ hit.prim = prim;
+ SET_PARAM(hit, st);
+ frag.time = time;
+ XD(setup_interface_fragment)(&frag, &vtx, &hit, fluid_side);
+ interf = scene_get_interface(scn, iprim);
+ ASSERT(interf);
+ *epsilon = interface_side_get_emissivity(interf, &frag);
+ #undef SET_PARAM
+ *hc = interface_get_convection_coef(interf, &frag);
+
+ return res;
+}
+
#if SDIS_SOLVE_DIMENSION == 3
static res_T
XD(ray_realisation)
@@ -1541,7 +1718,8 @@ XD(ray_realisation)
struct XD(temperature) T = XD(TEMPERATURE_NULL);
float dir[3];
res_T res = RES_OK;
- ASSERT(scn && position && direction && time>=0 && fp_to_meter>0 && weight);
+ ASSERT(scn && position && direction && time>=0 && fp_to_meter>0 && weight
+ && Tref >= 0);
ASSERT(medium && medium->type == SDIS_FLUID);
dX(set)(rwalk.vtx.P, position);
@@ -1610,6 +1788,11 @@ XD(solve_boundary)
SXD(scene_view_primitives_count(scn->sXd(view), &view_nprims));
FOR_EACH(i, 0, nprimitives) {
if(primitives[i] >= view_nprims) {
+ log_err(scn->dev,
+ "%s: invalid primitive identifier `%lu'. It must be in the [0 %lu] range.\n",
+ FUNC_NAME,
+ (unsigned long)primitives[i],
+ (unsigned long)scene_get_primitives_count(scn)-1);
res = RES_BAD_ARG;
goto error;
}
@@ -1617,9 +1800,9 @@ XD(solve_boundary)
/* Create the Star-XD shape of the boundary */
#if DIM == 2
- res = sXd(shape_create_line_segments)(scn->dev->sXd_dev, &shape);
+ res = s2d_shape_create_line_segments(scn->dev->sXd_dev, &shape);
#else
- res = sXd(shape_create_mesh)(scn->dev->sXd_dev, &shape);
+ res = s3d_shape_create_mesh(scn->dev->sXd_dev, &shape);
#endif
if(res != RES_OK) goto error;
@@ -1628,14 +1811,15 @@ XD(solve_boundary)
ctx.primitives = primitives;
ctx.view = scn->sXd(view);
vdata.usage = SXD_POSITION;
- vdata.type = DIM == 2 ? SXD_FLOAT2 : SXD_FLOAT3;
vdata.get = XD(boundary_get_position);
#if DIM == 2
- res = sXd(line_segments_setup_indexed_vertices)(shape, (unsigned)nprimitives,
- XD(boundary_get_indices), (unsigned)(nprimitives*DIM), &vdata, 1, &ctx);
+ vdata.type = S2D_FLOAT2;
+ res = s2d_line_segments_setup_indexed_vertices(shape, (unsigned)nprimitives,
+ boundary_get_indices_2d, (unsigned)(nprimitives*2), &vdata, 1, &ctx);
#else /* DIM == 3 */
- res = sXd(mesh_setup_indexed_vertices)(shape, (unsigned)nprimitives,
- XD(boundary_get_indices), (unsigned)(nprimitives*DIM), &vdata, 1, &ctx);
+ vdata.type = S3D_FLOAT3;
+ res = s3d_mesh_setup_indexed_vertices(shape, (unsigned)nprimitives,
+ boundary_get_indices_3d, (unsigned)(nprimitives*3), &vdata, 1, &ctx);
#endif
if(res != RES_OK) goto error;
@@ -1661,7 +1845,7 @@ XD(solve_boundary)
}
/* Create the estimator */
- res = estimator_create(scn->dev, &estimator);
+ res = estimator_create(scn->dev, SDIS_TEMPERATURE_ESTIMATOR, &estimator);
if(res != RES_OK) goto error;
omp_set_num_threads((int)scn->dev->nthreads);
@@ -1739,6 +1923,222 @@ error:
goto exit;
}
+static res_T
+XD(solve_boundary_flux)
+ (struct sdis_scene* scn,
+ const size_t nrealisations, /* #realisations */
+ const size_t primitives [], /* List of boundary primitives to handle */
+ const size_t nprimitives, /* #primitives */
+ const double time, /* Observation time */
+ const double fp_to_meter, /* Scale from floating point units to meters */
+ const double Tarad, /* In Kelvin */
+ const double Tref, /* In Kelvin */
+ struct sdis_estimator** out_estimator)
+{
+ struct XD(boundary_context) ctx = XD(BOUNDARY_CONTEXT_NULL);
+ struct sXd(vertex_data) vdata = SXD_VERTEX_DATA_NULL;
+ struct sXd(scene)* scene = NULL;
+ struct sXd(shape)* shape = NULL;
+ struct sXd(scene_view)* view = NULL;
+ struct sdis_estimator* estimator = NULL;
+ struct ssp_rng_proxy* rng_proxy = NULL;
+ struct ssp_rng** rngs = NULL;
+ double weight_t = 0, sqr_weight_t = 0;
+ double weight_fc = 0, sqr_weight_fc = 0;
+ double weight_fr = 0, sqr_weight_fr = 0;
+ double weight_f = 0, sqr_weight_f = 0;
+ size_t i;
+ size_t N = 0; /* #realisations that do not fail */
+ size_t view_nprims;
+ int64_t irealisation;
+ ATOMIC res = RES_OK;
+
+ if(!scn || !nrealisations || nrealisations > INT64_MAX || !primitives
+ || !nprimitives || time < 0 || fp_to_meter < 0 || Tref < 0
+ || !out_estimator) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ SXD(scene_view_primitives_count(scn->sXd(view), &view_nprims));
+ FOR_EACH(i, 0, nprimitives) {
+ if(primitives[i] >= view_nprims) {
+ log_err(scn->dev,
+ "%s: invalid primitive identifier `%lu'. It must be in the [0 %lu] range.\n",
+ FUNC_NAME,
+ (unsigned long)primitives[i],
+ (unsigned long)scene_get_primitives_count(scn)-1);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ }
+
+ /* Create the Star-XD shape of the boundary */
+#if DIM == 2
+ res = s2d_shape_create_line_segments(scn->dev->s2d, &shape);
+#else
+ res = s3d_shape_create_mesh(scn->dev->s3d, &shape);
+#endif
+ if(res != RES_OK) goto error;
+
+ /* Initialise the boundary shape with the triangles/segments of the
+ * submitted primitives */
+ ctx.primitives = primitives;
+ ctx.view = scn->sXd(view);
+ vdata.get = XD(boundary_get_position);
+#if DIM == 2
+ vdata.usage = S2D_POSITION;
+ vdata.type = S2D_FLOAT2;
+ res = s2d_line_segments_setup_indexed_vertices(shape, (unsigned)nprimitives,
+ boundary_get_indices_2d, (unsigned)(nprimitives*2), &vdata, 1, &ctx);
+#else /* DIM == 3 */
+ vdata.usage = S3D_POSITION;
+ vdata.type = S3D_FLOAT3;
+ res = s3d_mesh_setup_indexed_vertices(shape, (unsigned)nprimitives,
+ boundary_get_indices_3d, (unsigned)(nprimitives*3), &vdata, 1, &ctx);
+#endif
+ if(res != RES_OK) goto error;
+
+ /* Create and setup the boundary Star-XD scene */
+ res = sXd(scene_create)(scn->dev->sXd_dev, &scene);
+ if(res != RES_OK) goto error;
+ res = sXd(scene_attach_shape)(scene, shape);
+ if(res != RES_OK) goto error;
+ res = sXd(scene_view_create)(scene, SXD_SAMPLE, &view);
+ if(res != RES_OK) goto error;
+
+ /* Create the proxy RNG */
+ res = ssp_rng_proxy_create(scn->dev->allocator, &ssp_rng_mt19937_64,
+ scn->dev->nthreads, &rng_proxy);
+ if(res != RES_OK) goto error;
+
+ /* Create the per thread RNG */
+ rngs = MEM_CALLOC(scn->dev->allocator, scn->dev->nthreads, sizeof(*rngs));
+ if(!rngs) { res = RES_MEM_ERR; goto error; }
+ FOR_EACH(i, 0, scn->dev->nthreads) {
+ res = ssp_rng_proxy_create_rng(rng_proxy, i, rngs + i);
+ if(res != RES_OK) goto error;
+ }
+
+ /* Create the estimator */
+ res = estimator_create(scn->dev, SDIS_FLUX_ESTIMATOR, &estimator);
+ if(res != RES_OK) goto error;
+
+ omp_set_num_threads((int)scn->dev->nthreads);
+ #pragma omp parallel for schedule(static) reduction(+:weight_t,sqr_weight_t,\
+ weight_fc,sqr_weight_fc,weight_fr,sqr_weight_fr,weight_f,sqr_weight_f,N)
+ for(irealisation = 0; irealisation < (int64_t)nrealisations; ++irealisation) {
+ const int ithread = omp_get_thread_num();
+ struct sXd(primitive) prim;
+ struct ssp_rng* rng = rngs[ithread];
+ 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 };
+ double epsilon, hc, hr;
+ size_t iprim;
+ double uv[DIM - 1];
+ float st[DIM - 1];
+ res_T res_local = RES_OK;
+
+ if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occurred */
+
+ /* Sample a position onto the boundary */
+#if DIM == 2
+ res_local = s2d_scene_view_sample
+ (view,
+ ssp_rng_canonical_float(rng),
+ ssp_rng_canonical_float(rng),
+ &prim, st);
+ uv[0] = (double)st[0];
+#else
+ res_local = s3d_scene_view_sample
+ (view,
+ ssp_rng_canonical_float(rng),
+ ssp_rng_canonical_float(rng),
+ ssp_rng_canonical_float(rng),
+ &prim, st);
+ d2_set_f2(uv, st);
+#endif
+ if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); continue; }
+
+ /* Map from boundary scene to sdis scene */
+ ASSERT(prim.prim_id < nprimitives);
+ iprim = primitives[prim.prim_id];
+
+ interf = scene_get_interface(scn, (unsigned long)iprim);
+ fmd = interface_get_medium(interf, SDIS_FRONT);
+ bmd = interface_get_medium(interf, SDIS_BACK);
+ if(!fmd || !bmd
+ || (!(fmd->type == SDIS_FLUID && bmd->type == SDIS_SOLID)
+ && !(fmd->type == SDIS_SOLID && bmd->type == SDIS_FLUID)))
+ {
+ ATOMIC_SET(&res, RES_BAD_ARG);
+ continue;
+ }
+ solid_side = (fmd->type == SDIS_SOLID) ? SDIS_FRONT : SDIS_BACK;
+ fluid_side = (fmd->type == SDIS_FLUID) ? SDIS_FRONT : SDIS_BACK;
+
+ res_local = XD(interface_get_hc_epsilon)(&hc, &epsilon, scn,
+ (unsigned long)iprim, uv, time, fluid_side);
+ if(res_local != RES_OK) {
+ ATOMIC_SET(&res, res_local);
+ continue;
+ }
+ hr = 4.0 * BOLTZMANN_CONSTANT * Tref * Tref * Tref * epsilon;
+
+ /* Fluid, Radiative and Solid temperatures */
+ res_local = XD(probe_flux_realisation)(scn, rng, iprim, uv, time,
+ solid_side, fp_to_meter, Tarad, Tref, hr > 0, hc > 0, T_brf);
+ if(res_local != RES_OK) {
+ if(res_local != RES_BAD_OP) {
+ ATOMIC_SET(&res, res_local);
+ continue;
+ }
+ } else {
+ 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_total = w_conv + w_rad;
+ weight_t += Tboundary;
+ sqr_weight_t += Tboundary * Tboundary;
+ weight_fc += w_conv;
+ sqr_weight_fc += w_conv * w_conv;
+ weight_fr += w_rad;
+ sqr_weight_fr += w_rad * w_rad;
+ weight_f += w_total;
+ sqr_weight_f += w_total * w_total;
+ ++N;
+ }
+ }
+ if (res != RES_OK) goto error;
+
+ setup_estimator(estimator, nrealisations, N, weight_t, sqr_weight_t);
+ setup_estimator_flux(estimator, FLUX_CONVECTIVE__, weight_fc, sqr_weight_fc);
+ setup_estimator_flux(estimator, FLUX_RADIATIVE__, weight_fr, sqr_weight_fr);
+ setup_estimator_flux(estimator, FLUX_TOTAL__, weight_f, sqr_weight_f);
+
+exit:
+ if(scene) SXD(scene_ref_put(scene));
+ if(shape) SXD(shape_ref_put(shape));
+ 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(rngs) {
+ FOR_EACH(i, 0, scn->dev->nthreads) { if(rngs[i]) SSP(rng_ref_put(rngs[i])); }
+ MEM_RM(scn->dev->allocator, rngs);
+ }
+ return (res_T)res;
+error:
+ if(estimator) {
+ SDIS(estimator_ref_put(estimator));
+ estimator = NULL;
+ }
+ goto exit;
+}
+
#undef SDIS_SOLVE_DIMENSION
#undef DIM
#undef sXd
diff --git a/src/test_sdis_solve_boundary.c b/src/test_sdis_solve_boundary.c
@@ -311,6 +311,7 @@ main(int argc, char** argv)
uv[0] = 0.5;
iprim = 3;
+ CHK(SOLVE(square_scn, N, 4, uv, INF, F, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
CHK(SOLVE(square_scn, N, iprim, uv, INF, F, 1.0, 0, 0, &estimator) == RES_OK);
CHK(sdis_scene_get_boundary_position(square_scn, iprim, uv, pos) == RES_OK);
printf("Boundary temperature of the square at (%g %g) = ", SPLIT2(pos));
@@ -360,26 +361,27 @@ main(int argc, char** argv)
prims[0] = 4;
CHK(SOLVE(square_scn, N, prims, sides, 1, INF, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
- ref = (ref + Tb) / 2;
-
- /* Average temperature on the left/right side of the box */
+ /* Average temperature on the left+right sides of the box */
prims[0] = 2;
prims[1] = 3;
prims[2] = 6;
prims[3] = 7;
+
+ ref = (ref + Tb) / 2;
+
CHK(SOLVE(box_scn, N, prims, sides, 4, INF, 1.0, 0, 0, &estimator) == RES_OK);
- printf("Average temperature of the right/left side of the box = ");
+ printf("Average temperature of the left+right sides of the box = ");
check_estimator(estimator, N, ref);
CHK(sdis_estimator_ref_put(estimator) == RES_OK);
- /* Average temperature on the left/right side of the square */
+ /* Average temperature on the left+right sides of the square */
prims[0] = 1;
prims[1] = 3;
CHK(SOLVE(square_scn, N, prims, sides, 2, INF, 1.0, 0, 0, &estimator) == RES_OK);
- printf("Average temperature of the right/left side of the square = ");
+ printf("Average temperature of the left+right sides of the square = ");
check_estimator(estimator, N, ref);
CHK(sdis_estimator_ref_put(estimator) == RES_OK);
- #undef sdis_solve_boundary
+ #undef SOLVE
CHK(sdis_scene_ref_put(box_scn) == RES_OK);
CHK(sdis_scene_ref_put(square_scn) == RES_OK);
diff --git a/src/test_sdis_solve_boundary_flux.c b/src/test_sdis_solve_boundary_flux.c
@@ -0,0 +1,439 @@
+/* Copyright (C) 2016-2018 |Meso|Star> (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "sdis.h"
+#include "test_sdis_utils.h"
+
+#include <rsys/math.h>
+
+ /*
+ * The scene is composed of a solid cube/square whose temperature is unknown.
+ * The convection coefficient with the surrounding fluid is null excepted for
+ * the X faces whose value is 'H'. The Temperature T of the -X face is fixed
+ * to Tb. The ambiant radiative temperature is 0 excepted for the X faces
+ * whose value is 'Trad'.
+ * This test computes temperature and fluxes on the X faces and check that
+ * they are equal to:
+ *
+ * T(+X) = (H*Tf + Hrad*Trad + LAMBDA/A * Tb) / (H+Hrad+LAMBDA/A)
+ * with Hrad = 4 * BOLTZMANN_CONSTANT * Tref^3 * epsilon
+ * T(-X) = Tb
+ *
+ * CF = H * (T - Tf)
+ * RF = Hrad * (T - Trad)
+ * TF = CF + RF
+ *
+ * with Tf the temperature of the surrounding fluid, lambda the conductivity of
+ * the cube and A the size of the cube/square, i.e. 1.
+ *
+ * 3D
+ *
+ * ///////(1,1,1)
+ * +-------+
+ * /' /| _\ <-----
+ * -----> _\ +-------+ | / / H,Tf <----- Trad
+ * Trad -----> H,Tf / / Tb +.....|.+ \__/ <-----
+ * -----> \__/ |, |/
+ * +-------+
+ * (0,0,0)///////
+ *
+ *
+ * 2D
+ *
+ * ///////(1,1)
+ * +-------+
+ * -----> _\ | | _\ <-----
+ * Trad -----> H,Tf / / Tb | / / H,Tf <----- Trad
+ * -----> \__/ | | \__/ <-----
+ * +-------+
+ * (0,0)///////
+ */
+
+#define UNKNOWN_TEMPERATURE -1
+#define N 10000 /* #realisations */
+
+#define Tf 300.0
+#define Tb 0.0
+#define H 0.5
+#define Trad 300.0
+#define LAMBDA 0.1
+#define EPSILON 1.0
+
+#define Tref 300.0
+#define Hrad (4 * BOLTZMANN_CONSTANT * Tref * Tref * Tref * EPSILON)
+
+/*******************************************************************************
+ * Media
+ ******************************************************************************/
+struct fluid {
+ double temperature;
+};
+
+static double
+fluid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct fluid*)sdis_data_cget(data))->temperature;
+}
+
+
+static double
+solid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void) data;
+ CHK(vtx != NULL);
+ return 2.0;
+}
+
+static double
+solid_get_thermal_conductivity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void) data;
+ CHK(vtx != NULL);
+ return LAMBDA;
+}
+
+static double
+solid_get_volumic_mass
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void) data;
+ CHK(vtx != NULL);
+ return 25.0;
+}
+
+static double
+solid_get_delta
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void) data;
+ CHK(vtx != NULL);
+ return 1.0 / 20.0;
+}
+
+static double
+solid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void) data;
+ CHK(vtx != NULL);
+ return UNKNOWN_TEMPERATURE;
+}
+
+/*******************************************************************************
+ * Interfaces
+ ******************************************************************************/
+struct interf {
+ double temperature;
+ double emissivity;
+ double hc;
+};
+
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ const struct interf* interf = sdis_data_cget(data);
+ CHK(frag && data);
+ return interf->temperature;
+}
+
+static double
+interface_get_emissivity
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ const struct interf* interf = sdis_data_cget(data);
+ CHK(frag && data);
+ return interf->emissivity;
+}
+
+static double
+interface_get_convection_coef
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ const struct interf* interf = sdis_data_cget(data);
+ CHK(frag && data);
+ return interf->hc;
+}
+
+/*******************************************************************************
+ * Helper function
+ ******************************************************************************/
+static void
+check_estimator
+ (const struct sdis_estimator* estimator,
+ const size_t nrealisations, /* #realisations */
+ const double T,
+ const double CF,
+ const double RF,
+ const double TF)
+{
+ struct sdis_mc V = SDIS_MC_NULL;
+ enum sdis_estimator_type type;
+ size_t nreals;
+ size_t nfails;
+ CHK(estimator && nrealisations);
+
+ CHK(sdis_estimator_get_temperature(estimator, &V) == RES_OK);
+ CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
+ CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
+ printf("T = %g ~ %g +/- %g\n", T, V.E, V.SE);
+ CHK(eq_eps(V.E, T, 3 * (V.SE ? V.SE : FLT_EPSILON)));
+ CHK(sdis_estimator_get_type(estimator, &type) == RES_OK);
+ if(type == SDIS_FLUX_ESTIMATOR) {
+ CHK(sdis_estimator_get_convective_flux(estimator, &V) == RES_OK);
+ printf("Convective flux = %g ~ %g +/- %g\n", CF, V.E, V.SE);
+ CHK(eq_eps(V.E, CF, 3 * (V.SE ? V.SE : FLT_EPSILON)));
+ CHK(sdis_estimator_get_radiative_flux(estimator, &V) == RES_OK);
+ printf("Radiative flux = %g ~ %g +/- %g\n", RF, V.E, V.SE);
+ CHK(eq_eps(V.E, RF, 3 * (V.SE ? V.SE : FLT_EPSILON)));
+ CHK(sdis_estimator_get_total_flux(estimator, &V) == RES_OK);
+ printf("Total flux = %g ~ %g +/- %g\n", TF, V.E, V.SE);
+ CHK(eq_eps(V.E, TF, 3 * (V.SE ? V.SE : FLT_EPSILON)));
+ }
+ printf("#failures = %lu/%lu\n",
+ (unsigned long) nfails, (unsigned long) nrealisations);
+ CHK(nfails + nreals == nrealisations);
+ CHK(nfails < N / 1000);
+}
+
+/*******************************************************************************
+ * Test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct sdis_data* data = NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_medium* fluid = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_interface* interf_adiabatic = NULL;
+ struct sdis_interface* interf_Tb = NULL;
+ struct sdis_interface* interf_H = NULL;
+ struct sdis_scene* box_scn = NULL;
+ struct sdis_scene* square_scn = NULL;
+ struct sdis_estimator* estimator = NULL;
+ struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
+ struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
+ struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL;
+ struct sdis_interface* box_interfaces[12 /*#triangles*/];
+ struct sdis_interface* square_interfaces[4/*#segments*/];
+ struct interf* interf_props = NULL;
+ struct fluid* fluid_param;
+ enum sdis_estimator_type type;
+ double uv[2];
+ double pos[3];
+ double analyticT, analyticCF, analyticRF, analyticTF;
+ size_t prims[2];
+ size_t iprim;
+ (void)argc, (void)argv;
+
+ CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
+ CHK(sdis_device_create
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK);
+
+ /* Create the fluid medium */
+ CHK(sdis_data_create
+ (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data) == RES_OK);
+ fluid_param = sdis_data_get(data);
+ fluid_param->temperature = Tf;
+ fluid_shader.temperature = fluid_get_temperature;
+ CHK(sdis_fluid_create(dev, &fluid_shader, data, &fluid) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the solid_medium */
+ solid_shader.calorific_capacity = solid_get_calorific_capacity;
+ solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
+ solid_shader.volumic_mass = solid_get_volumic_mass;
+ solid_shader.delta_solid = solid_get_delta;
+ solid_shader.temperature = solid_get_temperature;
+ CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
+
+ /* Setup the interface shader */
+ interf_shader.convection_coef = interface_get_convection_coef;
+ interf_shader.front.temperature = interface_get_temperature;
+ interf_shader.front.specular_fraction = NULL;
+ interf_shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL;
+
+ /* Create the adiabatic interface */
+ CHK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data) == RES_OK);
+ interf_props = sdis_data_get(data);
+ interf_props->hc = 0;
+ interf_props->temperature = UNKNOWN_TEMPERATURE;
+ interf_props->emissivity = 0;
+ CHK(sdis_interface_create
+ (dev, solid, fluid, &interf_shader, data, &interf_adiabatic) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the Tb interface */
+ CHK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data) == RES_OK);
+ interf_props = sdis_data_get(data);
+ interf_props->hc = H;
+ interf_props->temperature = Tb;
+ interf_props->emissivity = EPSILON;
+ interf_shader.back.emissivity = interface_get_emissivity;
+ CHK(sdis_interface_create
+ (dev, solid, fluid, &interf_shader, data, &interf_Tb) == RES_OK);
+ interf_shader.back.emissivity = NULL;
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the H interface */
+ CHK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data) == RES_OK);
+ interf_props = sdis_data_get(data);
+ interf_props->hc = H;
+ interf_props->temperature = UNKNOWN_TEMPERATURE;
+ interf_props->emissivity = EPSILON;
+ interf_shader.back.emissivity = interface_get_emissivity;
+ CHK(sdis_interface_create
+ (dev, solid, fluid, &interf_shader, data, &interf_H) == RES_OK);
+ interf_shader.back.emissivity = NULL;
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Release the media */
+ CHK(sdis_medium_ref_put(solid) == RES_OK);
+ CHK(sdis_medium_ref_put(fluid) == RES_OK);
+
+ /* Map the interfaces to their box triangles */
+ box_interfaces[0] = box_interfaces[1] = interf_adiabatic; /* Front */
+ box_interfaces[2] = box_interfaces[3] = interf_Tb; /* Left */
+ box_interfaces[4] = box_interfaces[5] = interf_adiabatic; /* Back */
+ box_interfaces[6] = box_interfaces[7] = interf_H; /* Right */
+ box_interfaces[8] = box_interfaces[9] = interf_adiabatic; /* Top */
+ box_interfaces[10] = box_interfaces[11] = interf_adiabatic; /* Bottom */
+
+ /* Map the interfaces to their square segments */
+ square_interfaces[0] = interf_adiabatic; /* Bottom */
+ square_interfaces[1] = interf_Tb; /* Lef */
+ square_interfaces[2] = interf_adiabatic; /* Top */
+ square_interfaces[3] = interf_H; /* Right */
+
+ /* Create the box scene */
+ CHK(sdis_scene_create(dev, box_ntriangles, box_get_indices,
+ box_get_interface, box_nvertices, box_get_position, box_interfaces,
+ &box_scn) == RES_OK);
+
+ /* Create the square scene */
+ CHK(sdis_scene_2d_create(dev, square_nsegments, square_get_indices,
+ square_get_interface, square_nvertices, square_get_position,
+ square_interfaces, &square_scn) == RES_OK);
+
+ /* Release the interfaces */
+ CHK(sdis_interface_ref_put(interf_adiabatic) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_Tb) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_H) == RES_OK);
+
+ analyticT = (H*Tf + Hrad*Trad + LAMBDA * Tb) / (H + Hrad + LAMBDA);
+ analyticCF = H * (analyticT - Tf);
+ analyticRF = Hrad * (analyticT - Trad);
+ analyticTF = analyticCF + analyticRF;
+
+ #define SOLVE sdis_solve_probe_boundary_flux
+ uv[0] = 0.3;
+ uv[1] = 0.3;
+ iprim = 6;
+
+ CHK(SOLVE(NULL, N, iprim, uv, INF, 1.0, Trad, Tref, &estimator) == RES_BAD_ARG);
+ CHK(SOLVE(box_scn, 0, iprim, uv, INF, 1.0, Trad, Tref, &estimator) == RES_BAD_ARG);
+ CHK(SOLVE(box_scn, N, 12, uv, INF, 1.0, Trad, Tref, &estimator) == RES_BAD_ARG);
+ CHK(SOLVE(box_scn, N, iprim, NULL, INF, 1.0, Trad, Tref, &estimator) == RES_BAD_ARG);
+ CHK(SOLVE(box_scn, N, iprim, uv, -1, 1.0, Trad, Tref, &estimator) == RES_BAD_ARG);
+ CHK(SOLVE(box_scn, N, iprim, uv, INF, 1.0, Trad, Tref, NULL) == RES_BAD_ARG);
+
+ CHK(SOLVE(box_scn, N, iprim, uv, INF, 1.0, Trad, Tref, &estimator) == RES_OK);
+ CHK(sdis_estimator_get_type(estimator, &type) == RES_OK);
+ CHK(type == SDIS_FLUX_ESTIMATOR);
+
+ CHK(sdis_scene_get_boundary_position(box_scn, iprim, uv, pos) == RES_OK);
+ printf("Boundary values of the box at (%g %g %g) = ", SPLIT3(pos));
+ check_estimator(estimator, N, analyticT, analyticCF, analyticRF, analyticTF);
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+
+ uv[0] = 0.5;
+ iprim = 3;
+ CHK(SOLVE(square_scn, N, 4, uv, INF, 1.0, Trad, Tref, &estimator) == RES_BAD_ARG);
+ CHK(SOLVE(square_scn, N, iprim, uv, INF, 1.0, Trad, Tref, &estimator) == RES_OK);
+ CHK(sdis_scene_get_boundary_position(square_scn, iprim, uv, pos) == RES_OK);
+ printf("Boundary values of the square at (%g %g) = ", SPLIT2(pos));
+ check_estimator(estimator, N, analyticT, analyticCF, analyticRF, analyticTF);
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+
+ #undef F
+ #undef SOLVE
+
+ #define SOLVE sdis_solve_boundary_flux
+ prims[0] = 6;
+ prims[1] = 7;
+ CHK(SOLVE(NULL, N, prims, 2, INF, 1.0, Trad, Tref, &estimator) == RES_BAD_ARG);
+ CHK(SOLVE(box_scn, 0, prims, 2, INF, 1.0, Trad, Tref, &estimator) == RES_BAD_ARG);
+ CHK(SOLVE(box_scn, N, NULL, 2, INF, 1.0, Trad, Tref, &estimator) == RES_BAD_ARG);
+ CHK(SOLVE(box_scn, N, prims, 0, INF, 1.0, Trad, Tref, &estimator) == RES_BAD_ARG);
+ CHK(SOLVE(box_scn, N, prims, 2, -1, 1.0, Trad, Tref, &estimator) == RES_BAD_ARG);
+ CHK(SOLVE(box_scn, N, prims, 2, INF, 1.0, Trad, Tref, NULL) == RES_BAD_ARG);
+
+ /* Average temperature on the right side of the box */
+ CHK(SOLVE(box_scn, N, prims, 2, INF, 1.0, Trad, Tref, &estimator) == RES_OK);
+ printf("Average values of the right side of the box = ");
+ check_estimator(estimator, N, analyticT, analyticCF, analyticRF, analyticTF);
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+
+ /* Average temperature on the right side of the square */
+ prims[0] = 3;
+ CHK(SOLVE(square_scn, N, prims, 1, INF, 1.0, Trad, Tref, &estimator) == RES_OK);
+ printf("Average values of the right side of the square = ");
+ check_estimator(estimator, N, analyticT, analyticCF, analyticRF, analyticTF);
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+
+ /* Check out of bound prims */
+ prims[0] = 12;
+ CHK(SOLVE(box_scn, N, prims, 2, INF, 1.0, Trad, Tref, &estimator) == RES_BAD_ARG);
+ prims[0] = 4;
+ CHK(SOLVE(square_scn, N, prims, 1, INF, 1.0, Trad, Tref, &estimator) == RES_BAD_ARG);
+
+ /* Average temperature on the left side of the box */
+ prims[0] = 2;
+ prims[1] = 3;
+
+ analyticT = Tb;
+ analyticCF = H * (analyticT - Tf);
+ analyticRF = Hrad * (analyticT - Trad);
+ analyticTF = analyticCF + analyticRF;
+
+ CHK(SOLVE(box_scn, N, prims, 2, INF, 1.0, Trad, Tref, &estimator) == RES_OK);
+ printf("Average values of the left side of the box = ");
+ check_estimator(estimator, N, analyticT, analyticCF, analyticRF, analyticTF);
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+
+ /* Average temperature on the left/right side of the square */
+ prims[0] = 1;
+ CHK(SOLVE(square_scn, N, prims, 1, INF, 1.0, Trad, Tref, &estimator) == RES_OK);
+ printf("Average values of the left side of the square = ");
+ check_estimator(estimator, N, analyticT, analyticCF, analyticRF, analyticTF);
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+ #undef SOLVE
+
+ CHK(sdis_scene_ref_put(box_scn) == RES_OK);
+ CHK(sdis_scene_ref_put(square_scn) == RES_OK);
+ CHK(sdis_device_ref_put(dev) == RES_OK);
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
+
diff --git a/src/test_sdis_solve_probe.c b/src/test_sdis_solve_probe.c
@@ -174,6 +174,7 @@ main(int argc, char** argv)
{
struct mem_allocator allocator;
struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_mc F = SDIS_MC_NULL;
struct sdis_device* dev = NULL;
struct sdis_medium* solid = NULL;
struct sdis_medium* fluid = NULL;
@@ -188,6 +189,7 @@ main(int argc, char** argv)
struct fluid* fluid_param;
struct solid* solid_param;
struct interf* interface_param;
+ enum sdis_estimator_type type;
double pos[3];
double time;
double ref;
@@ -268,6 +270,24 @@ main(int argc, char** argv)
CHK(sdis_solve_probe(scn, N, pos, time, 1.0, 0, 0, NULL) == RES_BAD_ARG);
CHK(sdis_solve_probe(scn, N, pos, time, 1.0, 0, 0, &estimator) == RES_OK);
+ CHK(sdis_estimator_get_type(estimator, NULL) == RES_BAD_ARG);
+ CHK(sdis_estimator_get_type(NULL, &type) == RES_BAD_ARG);
+ CHK(sdis_estimator_get_type(estimator, &type) == RES_OK);
+ CHK(type == SDIS_TEMPERATURE_ESTIMATOR);
+
+ /* Fluxes aren't available after sdis_solve_probe */
+ CHK(sdis_estimator_get_convective_flux(estimator, NULL) == RES_BAD_ARG);
+ CHK(sdis_estimator_get_convective_flux(NULL, &F) == RES_BAD_ARG);
+ CHK(sdis_estimator_get_convective_flux(estimator, &F) == RES_BAD_ARG);
+
+ CHK(sdis_estimator_get_radiative_flux(estimator, NULL) == RES_BAD_ARG);
+ CHK(sdis_estimator_get_radiative_flux(NULL, &F) == RES_BAD_ARG);
+ CHK(sdis_estimator_get_radiative_flux(estimator, &F) == RES_BAD_ARG);
+
+ CHK(sdis_estimator_get_total_flux(estimator, NULL) == RES_BAD_ARG);
+ CHK(sdis_estimator_get_total_flux(NULL, &F) == RES_BAD_ARG);
+ CHK(sdis_estimator_get_total_flux(estimator, &F) == RES_BAD_ARG);
+
CHK(sdis_estimator_get_realisation_count(estimator, NULL) == RES_BAD_ARG);
CHK(sdis_estimator_get_realisation_count(NULL, &nreals) == RES_BAD_ARG);
CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
