stardis-solver

Solve coupled heat transfers
git clone git://git.meso-star.fr/stardis-solver.git
Log | Files | Refs | README | LICENSE

commit a46ab8f154b6596b815dd6aa8d8a3bfbc4de6b40
parent d0cb8e268b33426274697dd4b119b5a824b22f41
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Tue, 29 May 2018 12:34:34 +0200

Merge remote-tracking branch 'origin/test_volumic_power' into test_volumic_power

Diffstat:
Msrc/test_sdis_volumic_power3_2d.c | 253+++++++++++++++++++++++++++++++++++++++++++------------------------------------
1 file changed, 137 insertions(+), 116 deletions(-)

diff --git a/src/test_sdis_volumic_power3_2d.c b/src/test_sdis_volumic_power3_2d.c @@ -17,16 +17,40 @@ #include "test_sdis_utils.h" #include <rsys/math.h> -#define Pw 10000.0 -#define LAMBDA 10.0 -#define LAMBDA1 1.0 -#define LAMBDA2 LAMBDA1 -#define T1 373.15 -#define T2 273.15 -#define H1 5.0 -#define H2 10.0 -#define MDb 1.0 -#define N 400000 /* #realisations */ +#define Pw 10000.0 /* Volumic power */ +#define LAMBDA 10.0 /* Lambda of the middle slab */ +#define LAMBDA1 1.0 /* Lambda of the upper slab */ +#define LAMBDA2 LAMBDA1 /* Lambda of the lower slab */ +#define T1 373.15 /* Temperature of the upper fluid */ +#define T2 273.15 /* Temperature of the lower fluid */ +#define H1 5.0 /* Convection coef between the upper slab and the fluid */ +#define H2 10.0 /* Convection coef between the lower slab and the fluid */ +#define DELTA 0.01 /* Delta of the middle slab */ +#define DELTA1 0.02 /* Delta of the upper slab */ +#define DELTA2 0.07 /* Delta of the lower slab */ +#define MDb 1.0 /* Multiplier applied to delta to define delta boundary */ +#define L 0.2 /* Size of the middle slab */ +#define L1 0.4 /* Size of the upper slab */ +#define L2 1.4 /* Size of the lower slab */ + +#define N 10000 /* #realisations */ + +/* Analitically computed temperatures wrt the previous parameters.*/ +#define Tp1 648.6217 +#define Tp2 335.4141 +#define Ta 1199.5651 +#define Tb 1207.1122 + +/* Fixed temperatures */ +#define UNKNOWN_TEMPERATURE -1 +#define Tsolid1_fluid UNKNOWN_TEMPERATURE /*Tp1*/ +#define Tsolid2_fluid UNKNOWN_TEMPERATURE /*Tp2*/ +#define Tsolid_solid1 UNKNOWN_TEMPERATURE /*Ta*/ +#define Tsolid_solid2 UNKNOWN_TEMPERATURE /*Tb*/ + +/* Legacy deltas, 400K realisations: 924.093 ~ 928.405 +/- 0.847232 + * Deltas / 2, 400K realisations: 924.093 ~ 926.632 +/- 0.844205 */ +#define PROBE_POS 1.8 /* * The 2D scene is composed of 3 stacked solid slabs whose middle slab has a @@ -37,31 +61,31 @@ * _\ T1 * / / * \__/ - * ... -----H1------ ... + * ... -----H1------ ... Tp1 * LAMBDA1 * - * ... ------------- ... + * ... ------------- ... Ta * LAMBDA, Pw - * ... ------------- ... + * ... ------------- ... Tb * * LAMBDA2 * * - * ... -----H2------ ... + * ... -----H2------ ... Tp2 * _\ T2 * / / * \__/ */ static const double vertices[8/*#vertices*/*2/*#coords per vertex*/] = { - -100000.5, 0.0, - -100000.5, 1.4, - -100000.5, 1.6, - -100000.5, 2.0, - 100000.5, 2.0, - 100000.5, 1.6, - 100000.5, 1.4, - 100000.5, 0.0 + -100000.5, 0.0, /* 0 */ + -100000.5, 1.4, /* 1 */ + -100000.5, 1.6, /* 2 */ + -100000.5, 2.0, /* 3 */ + 100000.5, 2.0, /* 4 */ + 100000.5, 1.6, /* 5 */ + 100000.5, 1.4, /* 6 */ + 100000.5, 0.0 /* 7 */ }; static const size_t nvertices = sizeof(vertices)/sizeof(double[2]); @@ -232,23 +256,25 @@ main(int argc, char** argv) struct sdis_device* dev = NULL; struct sdis_data* data = NULL; struct sdis_medium* fluid = NULL; - struct sdis_medium* solid0 = NULL; + struct sdis_medium* solid = NULL; struct sdis_medium* solid1 = NULL; + struct sdis_medium* solid2 = NULL; struct sdis_scene* scn = NULL; struct sdis_estimator* estimator = NULL; struct sdis_fluid_shader fluid_shader = SDIS_FLUID_SHADER_NULL; struct sdis_solid_shader solid_shader = SDIS_SOLID_SHADER_NULL; struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; - struct sdis_interface* interf_adiabatic = NULL; - struct sdis_interface* interf_solid0_solid1_upp = NULL; - struct sdis_interface* interf_solid0_solid1_low = NULL; - struct sdis_interface* interf_solid0_upp = NULL; - struct sdis_interface* interf_solid0_low = NULL; + struct sdis_interface* interf_solid_adiabatic = NULL; struct sdis_interface* interf_solid1_adiabatic = NULL; + struct sdis_interface* interf_solid2_adiabatic = NULL; + struct sdis_interface* interf_solid_solid1 = NULL; + struct sdis_interface* interf_solid_solid2 = NULL; + struct sdis_interface* interf_solid1_fluid = NULL; + struct sdis_interface* interf_solid2_fluid = NULL; struct sdis_interface* interfaces[10/*#segment*/]; struct sdis_mc T = SDIS_MC_NULL; + double Tref; double pos[2]; - size_t i; (void)argc, (void)argv; CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK); @@ -276,112 +302,120 @@ main(int argc, char** argv) solid_shader.temperature = solid_get_temperature; solid_shader.volumic_power = solid_get_volumic_power; - /* Create the solid0 medium */ + /* Create the medium of the upper slab */ CHK(sdis_data_create (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data) == RES_OK); solid_param = sdis_data_get(data); solid_param->cp = 500000; solid_param->rho = 1000; solid_param->lambda = LAMBDA1; - solid_param->delta = 0.02; + solid_param->delta = DELTA1; solid_param->volumic_power = SDIS_VOLUMIC_POWER_NONE; solid_param->temperature = -1; - CHK(sdis_solid_create(dev, &solid_shader, data, &solid0) == RES_OK); + CHK(sdis_solid_create(dev, &solid_shader, data, &solid1) == RES_OK); CHK(sdis_data_ref_put(data) == RES_OK); - /* Create the solid1 medium */ + /* Create the medium of the lower slab */ + CHK(sdis_data_create + (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data) == RES_OK); + solid_param = sdis_data_get(data); + solid_param->cp = 500000; + solid_param->rho = 1000; + solid_param->lambda = LAMBDA2; + solid_param->delta = DELTA2; + solid_param->volumic_power = SDIS_VOLUMIC_POWER_NONE; + solid_param->temperature = -1; + CHK(sdis_solid_create(dev, &solid_shader, data, &solid2) == RES_OK); + CHK(sdis_data_ref_put(data) == RES_OK); + + /* Create the medium of the middle slab */ CHK(sdis_data_create (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data) == RES_OK); solid_param = sdis_data_get(data); solid_param->cp = 500000; solid_param->rho = 1000; solid_param->lambda = LAMBDA; - solid_param->delta = 0.01; + solid_param->delta = DELTA; solid_param->volumic_power = Pw; solid_param->temperature = -1; - CHK(sdis_solid_create(dev, &solid_shader, data, &solid1) == RES_OK); + CHK(sdis_solid_create(dev, &solid_shader, data, &solid) == RES_OK); CHK(sdis_data_ref_put(data) == RES_OK); interf_shader.front.temperature = interface_get_temperature; - /* Create the solid0/solid1 upper interface */ + /* Create the solid/solid1 interface */ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf), NULL, &data) == RES_OK); interf_param = sdis_data_get(data); - interf_param->temperature = -1/*1199.5651*/; - CHK(sdis_interface_create(dev, solid1, solid0, &interf_shader, - data, &interf_solid0_solid1_upp) == RES_OK); + interf_param->temperature = Tsolid_solid1; + CHK(sdis_interface_create(dev, solid, solid1, &interf_shader, + data, &interf_solid_solid1) == RES_OK); CHK(sdis_data_ref_put(data) == RES_OK); - /* Create the solid0/solid1 lower interface */ + /* Create the solid/solid2 interface */ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf), NULL, &data) == RES_OK); interf_param = sdis_data_get(data); - interf_param->temperature = -1/*1207.1122*/; - CHK(sdis_interface_create(dev, solid1, solid0, &interf_shader, - data, &interf_solid0_solid1_low) == RES_OK); + interf_param->temperature = Tsolid_solid2; + CHK(sdis_interface_create(dev, solid, solid2, &interf_shader, + data, &interf_solid_solid2) == RES_OK); CHK(sdis_data_ref_put(data) == RES_OK); /* Setup the interface shader */ interf_shader.convection_coef = interface_get_convection_coef; interf_shader.front.temperature = interface_get_temperature; - /* Create the adiabatic interface */ + /* Create the adiabatic interfaces */ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf), NULL, &data) == RES_OK); interf_param = sdis_data_get(data); interf_param->h = 0; interf_param->temperature = -1; - CHK(sdis_interface_create(dev, solid0, fluid, &interf_shader, data, - &interf_adiabatic) == RES_OK); - CHK(sdis_data_ref_put(data) == RES_OK); - - /* Create the solid0 fluid lower interface */ - CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf), - NULL, &data) == RES_OK); - interf_param = sdis_data_get(data); - interf_param->h = H2; - interf_param->temperature = 335.4141; - CHK(sdis_interface_create(dev, solid0, fluid, &interf_shader, data, - &interf_solid0_low) == RES_OK); + CHK(sdis_interface_create(dev, solid, fluid, &interf_shader, data, + &interf_solid_adiabatic) == RES_OK); + CHK(sdis_interface_create(dev, solid1, fluid, &interf_shader, data, + &interf_solid1_adiabatic) == RES_OK); + CHK(sdis_interface_create(dev, solid2, fluid, &interf_shader, data, + &interf_solid2_adiabatic) == RES_OK); CHK(sdis_data_ref_put(data) == RES_OK); - /* Create the solid0 upp interace */ + /* Create the solid1 fluid interace */ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf), NULL, &data) == RES_OK); interf_param = sdis_data_get(data); interf_param->h = H1; - interf_param->temperature = 648.6217; - CHK(sdis_interface_create(dev, solid0, fluid, &interf_shader, data, - &interf_solid0_upp) == RES_OK); + interf_param->temperature = Tsolid1_fluid; + CHK(sdis_interface_create(dev, solid1, fluid, &interf_shader, data, + &interf_solid1_fluid) == RES_OK); CHK(sdis_data_ref_put(data) == RES_OK); - /* Create the solid1 adiabatic interface */ + /* Create the solid2 fluid interface */ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf), NULL, &data) == RES_OK); interf_param = sdis_data_get(data); - interf_param->h = 0; - interf_param->temperature = -1; - CHK(sdis_interface_create(dev, solid1, fluid, &interf_shader, data, - &interf_solid1_adiabatic) == RES_OK); + interf_param->h = H2; + interf_param->temperature = Tsolid2_fluid; + CHK(sdis_interface_create(dev, solid2, fluid, &interf_shader, data, + &interf_solid2_fluid) == RES_OK); CHK(sdis_data_ref_put(data) == RES_OK); /* Release the media */ CHK(sdis_medium_ref_put(fluid) == RES_OK); - CHK(sdis_medium_ref_put(solid0) == RES_OK); + CHK(sdis_medium_ref_put(solid) == RES_OK); CHK(sdis_medium_ref_put(solid1) == RES_OK); + CHK(sdis_medium_ref_put(solid2) == RES_OK); /* Map the interfaces to their square segments */ - interfaces[0] = interf_adiabatic; - interfaces[1] = interf_solid1_adiabatic; - interfaces[2] = interf_adiabatic; - interfaces[3] = interf_solid0_upp; - interfaces[4] = interf_adiabatic; - interfaces[5] = interf_solid1_adiabatic; - interfaces[6] = interf_adiabatic; - interfaces[7] = interf_solid0_low; - interfaces[8] = interf_solid0_solid1_low; - interfaces[9] = interf_solid0_solid1_upp; + interfaces[0] = interf_solid2_adiabatic; + interfaces[1] = interf_solid_adiabatic; + interfaces[2] = interf_solid1_adiabatic; + interfaces[3] = interf_solid1_fluid; + interfaces[4] = interf_solid1_adiabatic; + interfaces[5] = interf_solid_adiabatic; + interfaces[6] = interf_solid2_adiabatic; + interfaces[7] = interf_solid2_fluid; + interfaces[8] = interf_solid_solid2; + interfaces[9] = interf_solid_solid1; #if 0 dump_segments(stdout, vertices, nvertices, indices, nsegments); @@ -393,49 +427,36 @@ main(int argc, char** argv) nvertices, get_position, interfaces, &scn) == RES_OK); /* Release the interfaces */ - CHK(sdis_interface_ref_put(interf_adiabatic) == RES_OK); - CHK(sdis_interface_ref_put(interf_solid0_upp) == RES_OK); - CHK(sdis_interface_ref_put(interf_solid0_low) == RES_OK); + CHK(sdis_interface_ref_put(interf_solid_adiabatic) == RES_OK); CHK(sdis_interface_ref_put(interf_solid1_adiabatic) == RES_OK); - CHK(sdis_interface_ref_put(interf_solid0_solid1_upp) == RES_OK); - CHK(sdis_interface_ref_put(interf_solid0_solid1_low) == RES_OK); - - FOR_EACH(i, 0, 8) { - const double l = 0.2; /* Size of the middle slab */ - const double l1 = 0.4; /* Size of the upper slab */ - const double l2 = 1.4; /* Size of the lower slab */ - double ta, tb; - double tp1, tp2; - double Tref; - - pos[0] = 0; - pos[1] = 0.7; /*1.85 - (double)i*0.2;*/ - - ta = 1199.5651; - tb = 1207.1122; - tp1 = 648.6217; - tp2 = 335.4141; - - if(pos[1] > 0 && pos[1] < l2) { /* Lower slab */ - Tref = tp2 + (tb - tp2) * pos[1] / l2; - } else if(pos[1] > l2 && pos[1] < l2 + l) { /* Middle slab */ - Tref = - (ta + tb) / 2 - + (ta - tb)/l * (pos[1] - (l2+l/2)) - + Pw * (l*l/4.0 - pow((pos[1] - (l2+l/2)), 2)) / (2*LAMBDA); - } else if(pos[1] > l2 + l && pos[1] < l2 + l1 + l) { - Tref = ta + (tp1 - ta) / l1 * (pos[1] - (l+l2)); - } else { - FATAL("Unreachable code.\n"); - } - - CHK(sdis_solve_probe(scn, N, pos, INF, 1.f, -1, 0, &estimator) == RES_OK); - CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK); - printf("Temperature at (%g %g) = %g ~ %g +/- %g\n", - SPLIT2(pos), Tref, T.E, T.SE); - CHK(sdis_estimator_ref_put(estimator) == RES_OK); + CHK(sdis_interface_ref_put(interf_solid2_adiabatic) == RES_OK); + CHK(sdis_interface_ref_put(interf_solid_solid1) == RES_OK); + CHK(sdis_interface_ref_put(interf_solid_solid2) == RES_OK); + CHK(sdis_interface_ref_put(interf_solid1_fluid) == RES_OK); + CHK(sdis_interface_ref_put(interf_solid2_fluid) == RES_OK); + + pos[0] = 0; + pos[1] = PROBE_POS; + + if(pos[1] > 0 && pos[1] < L2) { /* Lower slab */ + Tref = Tp2 + (Tb - Tp2) * pos[1] / L2; + } else if(pos[1] > L2 && pos[1] < L2 + L) { /* Middle slab */ + Tref = + (Ta + Tb) / 2 + + (Ta - Tb)/L * (pos[1] - (L2+L/2)) + + Pw * (L*L/4.0 - pow((pos[1] - (L2+L/2)), 2)) / (2*LAMBDA); + } else if(pos[1] > L2 + L && pos[1] < L2 + L1 + L) { + Tref = Ta + (Tp1 - Ta) / L1 * (pos[1] - (L+L2)); + } else { + FATAL("Unreachable code.\n"); } + CHK(sdis_solve_probe(scn, N, pos, INF, 1.f, -1, 0, &estimator) == RES_OK); + CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK); + printf("Temperature at (%g %g) = %g ~ %g +/- %g [%g, %g]\n", + SPLIT2(pos), Tref, T.E, T.SE, T.E-3*T.SE, T.E+3*T.SE); + CHK(sdis_estimator_ref_put(estimator) == RES_OK); + CHK(sdis_scene_ref_put(scn) == RES_OK); CHK(sdis_device_ref_put(dev) == RES_OK);