commit e1eb199ff4c694bcafcb9f26ba68818dad4dd5ab
parent 6faf87c3b3fa119309e9ab0b4d53b34cdf882ed3
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Fri, 6 Jul 2018 12:17:52 +0200
Add non uniform convection coefficients
Diffstat:
14 files changed, 532 insertions(+), 85 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -29,11 +29,11 @@ CMAKE_DEPENDENT_OPTION(ALL_TESTS
# Check dependencies
################################################################################
find_package(RCMake 0.3 REQUIRED)
-find_package(Star2D 0.1 REQUIRED)
+find_package(Star2D 0.1.1 REQUIRED)
find_package(Star3D 0.4 REQUIRED)
find_package(StarSP 0.7 REQUIRED)
-find_package(StarEnc 0.1 REQUIRED)
-find_package(StarEnc2D 0.1 REQUIRED)
+find_package(StarEnc 0.1.1 REQUIRED)
+find_package(StarEnc2D 0.1.1 REQUIRED)
find_package(RSys 0.6 REQUIRED)
find_package(OpenMP 1.2 REQUIRED)
@@ -141,6 +141,7 @@ if(NOT NO_TEST)
new_test(test_sdis_conducto_radiative)
new_test(test_sdis_conducto_radiative_2d)
new_test(test_sdis_convection)
+ new_test(test_sdis_convection_non_uniform)
new_test(test_sdis_data)
new_test(test_sdis_device)
new_test(test_sdis_flux)
diff --git a/src/sdis.h b/src/sdis.h
@@ -188,12 +188,15 @@ struct sdis_interface_shader {
/* May be NULL for solid/solid or if the convection coefficient is 0 onto
* the whole interface. */
sdis_interface_getter_T convection_coef; /* In W.K^-1.m^-2 */
+ /* Under no circumstance can convection_coef() return outside of
+ * [0 convection_coef_upper_bound] */
+ double convection_coef_upper_bound;
struct sdis_interface_side_shader front;
struct sdis_interface_side_shader back;
};
#define SDIS_INTERFACE_SHADER_NULL__ \
- {NULL, SDIS_INTERFACE_SIDE_SHADER_NULL__, SDIS_INTERFACE_SIDE_SHADER_NULL__}
+ {NULL, 0, SDIS_INTERFACE_SIDE_SHADER_NULL__, SDIS_INTERFACE_SIDE_SHADER_NULL__}
static const struct sdis_interface_shader SDIS_INTERFACE_SHADER_NULL =
SDIS_INTERFACE_SHADER_NULL__;
diff --git a/src/sdis_interface.c b/src/sdis_interface.c
@@ -54,6 +54,13 @@ check_interface_shader
"function of the interface shader should be NULL.\n", caller_name);
}
+ if(shader->convection_coef_upper_bound < 0) {
+ log_warn(dev,
+ "%s: Invalid upper bound for convection coefficient (%g).\n",
+ caller_name, shader->convection_coef_upper_bound);
+ if(type[0] == SDIS_FLUID || type[1] == SDIS_FLUID) return 0;
+ }
+
FOR_EACH(i, 0, 2) {
switch(type[i]) {
case SDIS_SOLID:
diff --git a/src/sdis_interface_c.h b/src/sdis_interface_c.h
@@ -70,6 +70,14 @@ interface_get_convection_coef
}
static INLINE double
+interface_get_convection_coef_upper_bound
+ (const struct sdis_interface* interf)
+{
+ ASSERT(interf);
+ return interf->shader.convection_coef_upper_bound;
+}
+
+static INLINE double
interface_side_get_temperature
(const struct sdis_interface* interf,
const struct sdis_interface_fragment* frag)
diff --git a/src/sdis_scene.c b/src/sdis_scene.c
@@ -103,6 +103,7 @@ scene_release(ref_T * ref)
darray_medium_release(&scn->media);
darray_prim_prop_release(&scn->prim_props);
htable_enclosure_release(&scn->enclosures);
+ htable_d_release(&scn->tmp_hc_ub);
if(scn->s2d_view) S2D(scene_view_ref_put(scn->s2d_view));
if(scn->s3d_view) S3D(scene_view_ref_put(scn->s3d_view));
MEM_RM(dev->allocator, scn);
diff --git a/src/sdis_scene_Xd.h b/src/sdis_scene_Xd.h
@@ -384,6 +384,8 @@ XD(setup_properties)
unsigned enclosures[2];
unsigned iprim_adjusted; /* Primitive id in user space */
unsigned id;
+ int i;
+ double* enc_upper_bound;
size_t ninterfaces;
#if DIM == 2
@@ -433,6 +435,17 @@ XD(setup_properties)
* front facing when their vertex are CCW ordered */
prim_prop->back_enclosure = enclosures[0];
prim_prop->front_enclosure = enclosures[1];
+
+ /* Build per-interface hc upper bounds in a tmp table */
+ FOR_EACH(i, 0, 2) {
+ enc_upper_bound = htable_d_find(&scn->tmp_hc_ub, enclosures+i);
+ double hc_ub = interface_get_convection_coef_upper_bound(itface);
+ if(!enc_upper_bound) {
+ res = htable_d_set(&scn->tmp_hc_ub, enclosures+i, &hc_ub);
+ } else {
+ *enc_upper_bound = MMAX(*enc_upper_bound, hc_ub);
+ }
+ }
}
exit:
@@ -510,6 +523,7 @@ XD(setup_enclosure_geometry)(struct sdis_scene* scn, struct sencXd(enclosure)* e
struct enclosure enc_dummy;
struct enclosure* enc_data;
float S, V;
+ double* p_ub;
unsigned iprim, nprims, nverts;
#if DIM == 2
struct senc2d_enclosure_header header;
@@ -581,6 +595,11 @@ XD(setup_enclosure_geometry)(struct sdis_scene* scn, struct sencXd(enclosure)* e
enc_data->S_over_V = S/absf(V);
#undef CALL
+ /* Set enclosure hc upper bound regardless of its media being a fluid */
+ p_ub = htable_d_find(&scn->tmp_hc_ub, &header.enclosure_id);
+ ASSERT(p_ub);
+ enc_data->hc_upper_bound = *p_ub;
+
/* Define the identifier of the enclosure primitives in the whole scene */
res = darray_uint_resize(&enc_data->local2global, nprims);
if(res != RES_OK) goto error;
@@ -639,6 +658,8 @@ XD(setup_enclosures)(struct sdis_scene* scn, struct sencXd(descriptor)* desc)
enc = NULL;
}
+ /* tmp table no more useful */
+ htable_d_purge(&scn->tmp_hc_ub);
exit:
if(enc) SENCXD(enclosure_ref_put(enc));
return res;
@@ -682,6 +703,7 @@ XD(scene_create)
darray_medium_init(dev->allocator, &scn->media);
darray_prim_prop_init(dev->allocator, &scn->prim_props);
htable_enclosure_init(dev->allocator, &scn->enclosures);
+ htable_d_init(dev->allocator, &scn->tmp_hc_ub);
res = XD(run_analyze)(scn, nprims, indices, interf, nverts, position, ctx, &desc);
if(res != RES_OK) {
diff --git a/src/sdis_scene_c.h b/src/sdis_scene_c.h
@@ -72,6 +72,7 @@ struct enclosure {
* whole scene */
struct darray_uint local2global;
+ double hc_upper_bound;
double S_over_V; /* in 3D = surface/volume; in 2D = perimeter/area */
};
@@ -83,6 +84,7 @@ enclosure_init(struct mem_allocator* allocator, struct enclosure* enc)
enc->s3d_view = NULL;
darray_uint_init(allocator, &enc->local2global);
enc->S_over_V = 0;
+ enc->hc_upper_bound = 0;
}
static INLINE void
@@ -105,6 +107,7 @@ enclosure_copy(struct enclosure* dst, const struct enclosure* src)
dst->s2d_view = src->s2d_view;
}
dst->S_over_V = src->S_over_V;
+ dst->hc_upper_bound = src->hc_upper_bound;
return darray_uint_copy(&dst->local2global, &src->local2global);
}
@@ -125,6 +128,7 @@ enclosure_copy_and_release(struct enclosure* dst, struct enclosure* src)
src->s2d_view = NULL;
}
dst->S_over_V = src->S_over_V;
+ dst->hc_upper_bound = src->hc_upper_bound;
return RES_OK;
}
@@ -156,6 +160,12 @@ enclosure_copy_and_release(struct enclosure* dst, struct enclosure* src)
#define HTABLE_DATA_FUNCTOR_COPY_AND_RELEASE enclosure_copy_and_release
#include <rsys/hash_table.h>
+/* Declare the hash table that maps an enclosure id to its data */
+#define HTABLE_NAME d
+#define HTABLE_KEY unsigned
+#define HTABLE_DATA double
+#include <rsys/hash_table.h>
+
struct sdis_scene {
struct darray_interf interfaces; /* List of interfaces own by the scene */
struct darray_medium media; /* List of media own by the scene */
@@ -163,6 +173,7 @@ struct sdis_scene {
struct s2d_scene_view* s2d_view;
struct s3d_scene_view* s3d_view;
+ struct htable_d tmp_hc_ub; /* Map an enclosure id to its hc upper bound */
struct htable_enclosure enclosures; /* Map an enclosure id to its data */
double ambient_radiative_temperature; /* In Kelvin */
diff --git a/src/sdis_solve_Xd.h b/src/sdis_solve_Xd.h
@@ -443,6 +443,7 @@ XD(fluid_temperature)
double mu;
double tau;
double tmp;
+ double r;
#if DIM == 2
float st;
#else
@@ -453,7 +454,7 @@ XD(fluid_temperature)
ASSERT(rwalk->mdm->type == SDIS_FLUID);
tmp = fluid_get_temperature(rwalk->mdm, &rwalk->vtx);
- if(tmp >= 0) {
+ if(tmp >= 0) { /* T is known. */
T->value += tmp;
T->done = 1;
return RES_OK;
@@ -474,16 +475,13 @@ XD(fluid_temperature)
if(SXD_HIT_NONE(&rwalk->hit)) {
log_err(scn->dev,
-"%s: the position %g %g %g lise in the surrounding fluid whose temperature must \n"
+"%s: the position %g %g %g lies in the surrounding fluid whose temperature must \n"
"be known.\n",
FUNC_NAME, SPLIT3(rwalk->vtx.P));
return RES_BAD_OP;
}
}
- /* Setup the fragment of the interface */
- XD(setup_interface_fragment)(&frag, &rwalk->vtx, &rwalk->hit, rwalk->hit_side);
-
/* Fetch the current interface and its associated enclosures */
interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
scene_get_enclosure_ids(scn, rwalk->hit.prim.prim_id, enc_ids);
@@ -492,9 +490,11 @@ XD(fluid_temperature)
ASSERT(interf->medium_front != interf->medium_back);
if(rwalk->mdm == interf->medium_front) {
enc_id = enc_ids[0];
+ ASSERT(rwalk->hit_side == SDIS_FRONT);
} else {
ASSERT(rwalk->mdm == interf->medium_back);
enc_id = enc_ids[1];
+ ASSERT(rwalk->hit_side == SDIS_BACK);
}
/* Fetch the enclosure data */
@@ -502,63 +502,113 @@ XD(fluid_temperature)
if(!enc) {
log_err(scn->dev,
"%s: invalid enclosure. The position %g %g %g may lie in the surrounding fluid.\n",
- FUNC_NAME, SPLIT3(rwalk->vtx.P));
+ FUNC_NAME, SPLIT3(rwalk->vtx.P));
return RES_BAD_OP;
}
- /* Fetch the physical properties */
- hc = interface_get_convection_coef(interf, &frag);
- cp = fluid_get_calorific_capacity(rwalk->mdm, &rwalk->vtx);
- rho = fluid_get_volumic_mass(rwalk->mdm, &rwalk->vtx);
-
- /* Sample the time.
- * FIXME we assume that hc is constant for the whole enclosure */
- mu = hc / (rho * cp) * enc->S_over_V;
- tau = ssp_ran_exp(rng, mu);
- rwalk->vtx.time = MMAX(rwalk->vtx.time - tau, 0);
-
- /* Check the initial condition */
- tmp = fluid_get_temperature(rwalk->mdm, &rwalk->vtx);
- if(tmp >= 0) {
- T->value += tmp;
- T->done = 1;
- return RES_OK;
- }
+ /* The hc upper bound can be 0 is h is uniformly 0.
+ * In that case the result is the initial condition. */
+ if(enc->hc_upper_bound == 0) {
+ /* Cannot be in the fluid without starting there. */
+ ASSERT(SXD_HIT_NONE(&rwalk->hit));
+ rwalk->vtx.time = 0;
+ tmp = fluid_get_temperature(rwalk->mdm, &rwalk->vtx);
+ if(tmp >= 0) {
+ T->value += tmp;
+ T->done = 1;
+ return RES_OK;
+ }
- /* The initial condition should be reached */
- if(rwalk->vtx.time <=0) {
+ /* At t=0, the initial condition should have been reached. */
log_err(scn->dev,
- "%s: undefined initial condition. "
- "The time is null but the temperature remains unknown.\n",
+"%s: undefined initial condition. "
+"Time is 0 but the temperature remains unknown.\n",
FUNC_NAME);
return RES_BAD_OP;
}
- /* Uniformly sample the enclosure */
+ /* A trick to force first r test result. */
+ r = 1;
+
+ /* Sample time until intial condition is reached
+ * or a true convection occurs. */
+ while(1) {
+ /* Setup the fragment of the interface. */
+ XD(setup_interface_fragment)(&frag, &rwalk->vtx, &rwalk->hit, rwalk->hit_side);
+
+ /* Fetch hc. */
+ hc = interface_get_convection_coef(interf, &frag);
+ if(hc > enc->hc_upper_bound) {
+ log_err(scn->dev,
+ "%s: hc (%g) exceeds its provided upper bound (%g) at %g %g %.\n",
+ FUNC_NAME, hc, enc->hc_upper_bound, SPLIT3(rwalk->vtx.P));
+ return RES_BAD_OP;
+ }
+
+ if(r < hc / enc->hc_upper_bound) {
+ /* True convection. Always true if hc == bound. */
+ break;
+ }
+
+ /* Fetch other physical properties. */
+ cp = fluid_get_calorific_capacity(rwalk->mdm, &rwalk->vtx);
+ rho = fluid_get_volumic_mass(rwalk->mdm, &rwalk->vtx);
+
+ /* Sample the time using the upper bound. */
+ mu = enc->hc_upper_bound / (rho * cp) * enc->S_over_V;
+ tau = ssp_ran_exp(rng, mu);
+ rwalk->vtx.time = MMAX(rwalk->vtx.time - tau, 0);
+
+ /* Check the initial condition. */
+ tmp = fluid_get_temperature(rwalk->mdm, &rwalk->vtx);
+ if(tmp >= 0) {
+ T->value += tmp;
+ T->done = 1;
+ return RES_OK;
+ }
+
+ if(rwalk->vtx.time <= 0) {
+ /* The initial condition should have been reached. */
+ log_err(scn->dev,
+"%s: undefined initial condition. "
+"Time is 0 but the temperature remains unknown.\n",
+ FUNC_NAME);
+ return RES_BAD_OP;
+ }
+
+ /* Uniformly sample the enclosure. */
#if DIM == 2
- SXD(scene_view_sample
+ SXD(scene_view_sample
(enc->sXd(view),
- ssp_rng_canonical_float(rng),
- ssp_rng_canonical_float(rng),
- &rwalk->hit.prim,
- &rwalk->hit.u));
- st = rwalk->hit.u;
+ ssp_rng_canonical_float(rng),
+ ssp_rng_canonical_float(rng),
+ &rwalk->hit.prim,
+ &rwalk->hit.u));
+ st = rwalk->hit.u;
#else
- SXD(scene_view_sample
+ SXD(scene_view_sample
(enc->sXd(view),
- ssp_rng_canonical_float(rng),
- ssp_rng_canonical_float(rng),
- ssp_rng_canonical_float(rng),
- &rwalk->hit.prim,
- rwalk->hit.uv));
- f2_set(st, rwalk->hit.uv);
+ ssp_rng_canonical_float(rng),
+ ssp_rng_canonical_float(rng),
+ ssp_rng_canonical_float(rng),
+ &rwalk->hit.prim,
+ rwalk->hit.uv));
+ f2_set(st, rwalk->hit.uv);
#endif
- rwalk->hit.distance = 0;
- SXD(primitive_get_attrib(&rwalk->hit.prim, SXD_POSITION, st, &attr_P));
- SXD(primitive_get_attrib(&rwalk->hit.prim, SXD_GEOMETRY_NORMAL, st, &attr_N));
- dX_set_fX(rwalk->vtx.P, attr_P.value);
- fX(set)(rwalk->hit.normal, attr_N.value);
+ SXD(primitive_get_attrib(&rwalk->hit.prim, SXD_POSITION, st, &attr_P));
+ SXD(primitive_get_attrib(&rwalk->hit.prim, SXD_GEOMETRY_NORMAL, st, &attr_N));
+ dX_set_fX(rwalk->vtx.P, attr_P.value);
+ fX(set)(rwalk->hit.normal, attr_N.value);
+
+ /* Fetch the interface of the sampled point. */
+ interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
+
+ /* Renew r for next loop. */
+ r = ssp_rng_canonical_float(rng);
+ }
+
+ rwalk->hit.distance = 0;
T->func = XD(boundary_temperature);
rwalk->mdm = NULL; /* The random walk is at an interface between 2 media */
return RES_OK;
diff --git a/src/test_sdis_conducto_radiative.c b/src/test_sdis_conducto_radiative.c
@@ -239,6 +239,7 @@ create_interface
shader.back.emissivity = interface_get_emissivity;
shader.back.specular_fraction = interface_get_specular_fraction;
}
+ shader.convection_coef_upper_bound = MMAX(0, interf->convection_coef);
CHK(sdis_data_create(dev, sizeof(struct interfac), ALIGNOF(struct interfac),
NULL, &data) == RES_OK);
diff --git a/src/test_sdis_conducto_radiative_2d.c b/src/test_sdis_conducto_radiative_2d.c
@@ -255,6 +255,8 @@ create_interface
shader.back.emissivity = interface_get_emissivity;
shader.back.specular_fraction = interface_get_specular_fraction;
}
+ shader.convection_coef_upper_bound = MMAX(0, interf->convection_coef);
+
CHK(sdis_data_create(dev, sizeof(struct interfac), ALIGNOF(struct interfac),
NULL, &data) == RES_OK);
*((struct interfac*)sdis_data_get(data)) = *interf;
diff --git a/src/test_sdis_convection.c b/src/test_sdis_convection.c
@@ -185,9 +185,9 @@ main(int argc, char** argv)
double ref;
double Tinf;
double nu;
- double time;
size_t nreals;
size_t nfails;
+ int i;
(void)argc, (void)argv;
CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
@@ -208,6 +208,7 @@ main(int argc, char** argv)
interf_shader.front.temperature = interface_get_temperature;
interf_shader.front.emissivity = interface_get_emissivity;
interf_shader.front.specular_fraction = interface_get_specular_fraction;
+ interf_shader.convection_coef_upper_bound = H;
/* Create the interfaces */
interf_T0 = create_interface(dev, fluid, solid, &interf_shader, T0);
@@ -255,39 +256,47 @@ main(int argc, char** argv)
d3_splat(pos, 0.25);
- nu = (6*H)/(RHO*CP);
- time = 1.0/nu;
- Tinf = (H*(T0+T1+T2+T3+T4+T5))/(6*H);
- ref = Tf_0 * exp(-nu*time) + Tinf*(1-exp(-nu*time));
-
- /* Solve in 3D */
- CHK(sdis_solve_probe(box_scn, N, pos, time, 1.0, 0, 0, &estimator) == RES_OK);
- CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
- CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
- CHK(nfails + nreals == N);
- CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
- CHK(sdis_estimator_ref_put(estimator) == RES_OK);
- printf("Temperature of the box at (%g %g %g) = %g ~ %g +/- %g\n",
- SPLIT3(pos), ref, T.E, T.SE);
- printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
- CHK(eq_eps(T.E, ref, T.SE*2));
-
- nu = (4*H)/(RHO*CP);
- time = 1.0/nu;
- Tinf = (H*(T0+T1+T2+T3))/(4*H);
- ref = Tf_0 * exp(-nu*time) + Tinf*(1-exp(-nu*time));
-
- /* Solve in 2D */
- CHK(sdis_solve_probe(square_scn, N, pos, time, 1.0, 0, 0, &estimator) == RES_OK);
- CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
- CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
- CHK(nfails + nreals == N);
- CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
- CHK(sdis_estimator_ref_put(estimator) == RES_OK);
- printf("Temperature of the square at (%g %g) = %g ~ %g +/- %g\n",
- SPLIT2(pos), ref, T.E, T.SE);
- printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
- CHK(eq_eps(T.E, ref, T.SE*2.0));
+ /* Test in 3D for various time values. */
+ nu = (6 * H) / (RHO*CP);
+ Tinf = (H*(T0 + T1 + T2 + T3 + T4 + T5)) / (6 * H);
+ printf("Temperature of the box at (%g %g %g)\n", SPLIT3(pos));
+ FOR_EACH(i, 0, 5) {
+ double time = i ? (double) i / nu : INF;
+ ref = Tf_0 * exp(-nu * time) + Tinf * (1 - exp(-nu * time));
+
+ /* Solve in 3D */
+ CHK(sdis_solve_probe(box_scn, N, pos, time, 1.0, 0, 0, &estimator) == RES_OK);
+ CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
+ CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
+ CHK(nfails + nreals == N);
+ CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+ printf(" t=%g : %g ~ %g +/- %g\n", time, ref, T.E, T.SE);
+ if(nfails)
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
+ CHK(eq_eps(T.E, ref, T.SE * 3));
+ }
+
+ /* Test in 2D for various time values. */
+ nu = (4 * H) / (RHO*CP);
+ Tinf = (H * (T0 + T1 + T2 + T3)) / (4 * H);
+ printf("Temperature of the square at (%g %g)\n", SPLIT2(pos));
+ FOR_EACH(i, 0, 5) {
+ double time = i ? (double) i / nu : INF;
+ ref = Tf_0 * exp(-nu * time) + Tinf * (1 - exp(-nu * time));
+
+ /* Solve in 2D */
+ CHK(sdis_solve_probe(square_scn, N, pos, time, 1.0, 0, 0, &estimator) == RES_OK);
+ CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
+ CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
+ CHK(nfails + nreals == N);
+ CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+ printf(" t=%g : %g ~ %g +/- %g\n", time, ref, T.E, T.SE);
+ if(nfails)
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
+ CHK(eq_eps(T.E, ref, T.SE * 3));
+ }
CHK(sdis_scene_ref_put(box_scn) == RES_OK);
CHK(sdis_scene_ref_put(square_scn) == RES_OK);
diff --git a/src/test_sdis_convection_non_uniform.c b/src/test_sdis_convection_non_uniform.c
@@ -0,0 +1,325 @@
+/* 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/double3.h>
+#include <rsys/math.h>
+
+/*
+ * The scene is composed of an unit fluid cube/square whose temperature is
+ * unknown. The convection coefficient with the surrounding solid is H
+ * everywhere the temperature of the -/+X, -/+Y and -/+Z faces are fixed to T0
+ * and T1, T2, T3, T4 and T5, respectively. This test computes the temperature
+ * of the fluid Tf at an observation time t. This temperature is equal to:
+ *
+ * Tf(t) = Tf(0) * e^(-nu*t) + Tinf*(1-e^(-nu*t))
+ *
+ * nu = (Sum_{i=0..5}(H*Si)) / (RHO*CP*V)
+ * Tinf = (Sum_{i=0..5}(H*Si*Ti) / (Sum_{i=0..5}(H*Si));
+ *
+ * with Si surface of the faces (i.e. one), V the volume of the cube (i.e.
+ * one), RHO the volumic mass of the fluid and CP its calorific capacity.
+ *
+ * 3D 2D
+ *
+ * (1,1,1) (1,1)
+ * +---------+ +-----T3----+
+ * /' T3 /|T4 | |
+ * +---------+ | | H _\ |
+ * | ' H _\ |T1 T0 / / T1
+ * |T0 / / | | | \__/ |
+ * | +..\__/.|.+ | |
+ * T5|, T2 |/ +-----------+
+ * +---------+ (0,0)
+ * (0,0,0)
+ */
+
+#define UNKNOWN_TEMPERATURE -1
+#define N 100000 /* #realisations */
+
+#define Tf_0 280.0
+
+#define T0 300.0
+#define T1 310.0
+#define T2 320.0
+#define T3 330.0
+#define T4 340.0
+#define T5 350.0
+
+#define HC0 100.0
+#define HC1 30.0
+#define HC2 3020.0
+#define HC3 7300.0
+#define HC4 3400.0
+#define HC5 50.0
+
+#define RHO 25.0
+#define CP 2.0
+
+/*******************************************************************************
+ * Media
+ ******************************************************************************/
+static double
+fluid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return vtx->time <= 0 ? Tf_0 : UNKNOWN_TEMPERATURE;
+}
+
+static double
+fluid_get_volumic_mass
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return RHO;
+}
+
+static double
+fluid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return CP;
+}
+
+/*******************************************************************************
+ * Interface
+ ******************************************************************************/
+struct interf {
+ double temperature;
+ 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_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;
+}
+
+static double
+interface_get_emissivity
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag && data);
+ return 0;
+}
+
+static double
+interface_get_specular_fraction
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag && data);
+ return 0;
+}
+
+static struct sdis_interface*
+create_interface
+ (struct sdis_device* dev,
+ struct sdis_medium* front,
+ struct sdis_medium* back,
+ const struct sdis_interface_shader* interf_shader,
+ const double temperature,
+ const double hc)
+{
+ struct sdis_data* data;
+ struct sdis_interface* interf;
+ struct interf* interf_props;
+
+ CHK(sdis_data_create
+ (dev, sizeof(struct interf), ALIGNOF(struct interf), NULL, &data) == RES_OK);
+ interf_props = sdis_data_get(data);
+ interf_props->temperature = temperature;
+ interf_props->hc = hc;
+ CHK(sdis_interface_create
+ (dev, front, back, interf_shader, data, &interf) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+ return interf;
+}
+
+/*******************************************************************************
+ * Test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_medium* fluid = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_interface* interf_T0 = NULL;
+ struct sdis_interface* interf_T1 = NULL;
+ struct sdis_interface* interf_T2 = NULL;
+ struct sdis_interface* interf_T3 = NULL;
+ struct sdis_interface* interf_T4 = NULL;
+ struct sdis_interface* interf_T5 = 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 = DUMMY_INTERFACE_SHADER;
+ struct sdis_interface* box_interfaces[12/*#triangles*/];
+ struct sdis_interface* square_interfaces[4/*#segments*/];
+ double pos[3];
+ double ref;
+ double Tinf;
+ double nu;
+ size_t nreals;
+ size_t nfails;
+ int i;
+ (void)argc, (void)argv;
+
+ CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
+ CHK(sdis_device_create
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 0, &dev) == RES_OK);
+
+ /* Create the fluid medium */
+ fluid_shader.temperature = fluid_get_temperature;
+ fluid_shader.calorific_capacity = fluid_get_calorific_capacity;
+ fluid_shader.volumic_mass = fluid_get_volumic_mass;
+ CHK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid) == RES_OK);
+
+ /* Create the solid_medium */
+ 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.emissivity = interface_get_emissivity;
+ interf_shader.front.specular_fraction = interface_get_specular_fraction;
+
+ /* Create the interfaces */
+ interf_shader.convection_coef_upper_bound = HC0;
+ interf_T0 = create_interface(dev, fluid, solid, &interf_shader, T0, HC0);
+ interf_shader.convection_coef_upper_bound = HC1;
+ interf_T1 = create_interface(dev, fluid, solid, &interf_shader, T1, HC1);
+ interf_shader.convection_coef_upper_bound = HC2;
+ interf_T2 = create_interface(dev, fluid, solid, &interf_shader, T2, HC2);
+ interf_shader.convection_coef_upper_bound = HC3;
+ interf_T3 = create_interface(dev, fluid, solid, &interf_shader, T3, HC3);
+ interf_shader.convection_coef_upper_bound = HC4;
+ interf_T4 = create_interface(dev, fluid, solid, &interf_shader, T4, HC4);
+ interf_shader.convection_coef_upper_bound = HC5;
+ interf_T5 = create_interface(dev, fluid, solid, &interf_shader, T5, HC5);
+
+ /* 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_T5; /* Front */
+ box_interfaces[2] = box_interfaces[3] = interf_T0; /* Left */
+ box_interfaces[4] = box_interfaces[5] = interf_T4; /* Back */
+ box_interfaces[6] = box_interfaces[7] = interf_T1; /* Right */
+ box_interfaces[8] = box_interfaces[9] = interf_T3; /* Top */
+ box_interfaces[10]= box_interfaces[11]= interf_T2; /* Bottom */
+
+ /* Map the interfaces to their square segments */
+ square_interfaces[0] = interf_T2; /* Bottom */
+ square_interfaces[1] = interf_T0; /* Left */
+ square_interfaces[2] = interf_T3; /* Top */
+ square_interfaces[3] = interf_T1; /* 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_T0) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_T1) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_T2) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_T3) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_T4) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_T5) == RES_OK);
+
+ d3_splat(pos, 0.25);
+
+ /* Test in 3D for various time values. */
+ nu = (HC0 + HC1 + HC2 + HC3 + HC4 + HC5) / (RHO * CP);
+ Tinf = (HC0 * T0 + HC1 * T1 + HC2 * T2 + HC3 * T3 + HC4 * T4 + HC5 * T5)
+ / (HC0 + HC1 + HC2 + HC3 + HC4 + HC5);
+ printf("Temperature of the box at (%g %g %g)\n", SPLIT3(pos));
+ FOR_EACH(i, 0, 5) {
+ double time = i ? (double) i / nu : INF;
+ ref = Tf_0 * exp(-nu * time) + Tinf * (1 - exp(-nu * time));
+
+ /* Solve in 3D */
+ CHK(sdis_solve_probe(box_scn, N, pos, time, 1.0, 0, 0, &estimator) == RES_OK);
+ CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
+ CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
+ CHK(nfails + nreals == N);
+ CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+ printf(" t=%g : %g ~ %g +/- %g\n", time, ref, T.E, T.SE);
+ if(nfails)
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails,(unsigned long)N);
+ CHK(eq_eps(T.E, ref, T.SE * 3));
+ }
+
+ /* Test in 2D for various time values. */
+ nu = (HC0 + HC1 + HC2 + HC3) / (RHO * CP);
+ Tinf = (HC0 * T0 + HC1 * T1 + HC2 * T2 + HC3 * T3) / (HC0 + HC1 + HC2 + HC3);
+ printf("Temperature of the square at (%g %g)\n", SPLIT2(pos));
+ FOR_EACH(i, 0, 5) {
+ double time = i ? (double) i / nu : INF;
+ ref = Tf_0 * exp(-nu * time) + Tinf * (1 - exp(-nu * time));
+
+ CHK(sdis_solve_probe(square_scn, N, pos, time, 1.0, 0, 0, &estimator) == RES_OK);
+ CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
+ CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
+ CHK(nfails + nreals == N);
+ CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+ printf(" t=%g : %g ~ %g +/- %g\n", time, ref, T.E, T.SE);
+ if (nfails)
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails,(unsigned long)N);
+ CHK(eq_eps(T.E, ref, T.SE * 3));
+ }
+
+ 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_interface.c b/src/test_sdis_interface.c
@@ -103,6 +103,12 @@ main(int argc, char** argv)
shader.front.specular_fraction = dummy_interface_getter;
CHK(CREATE(dev, solid, fluid, &shader, NULL, &interf) == RES_OK); /* Warning */
CHK(sdis_interface_ref_put(interf) == RES_OK);
+ shader.front.specular_fraction = NULL;
+ shader.convection_coef_upper_bound = -1;
+ CHK(CREATE(dev, solid, solid, &shader, NULL, &interf) == RES_OK); /* Warning */
+ CHK(sdis_interface_ref_put(interf) == RES_OK);
+ CHK(CREATE(dev, solid, fluid, &shader, NULL, &interf) == RES_BAD_ARG);
+ shader.convection_coef_upper_bound = 0;
#undef CREATE
CHK(sdis_device_ref_put(dev) == RES_OK);
diff --git a/src/test_sdis_utils.h b/src/test_sdis_utils.h
@@ -183,6 +183,7 @@ static const struct sdis_fluid_shader DUMMY_FLUID_SHADER = {
}
static const struct sdis_interface_shader DUMMY_INTERFACE_SHADER = {
dummy_interface_getter,
+ 0,
DUMMY_INTERFACE_SIDE_SHADER__,
DUMMY_INTERFACE_SIDE_SHADER__
};
