stardis-solver

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

commit 02661cd5733476c3866759359b2d058e2b875e15
parent 214747508884aaae7069c7861c4553e4bf693216
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Tue,  5 Jun 2018 18:23:04 +0200

Add some results to the volumic power tests

Diffstat:
Msrc/test_sdis_solve_probe_boundary.c | 5+++--
Msrc/test_sdis_volumic_power.c | 3++-
Msrc/test_sdis_volumic_power2.c | 12+++++++++++-
Msrc/test_sdis_volumic_power2_2d.c | 30++++++++++++++++++++++++++++--
Msrc/test_sdis_volumic_power4_2d.c | 107++++++++++++++++++++++++++++++++++++++++++++++++++++---------------------------
5 files changed, 114 insertions(+), 43 deletions(-)

diff --git a/src/test_sdis_solve_probe_boundary.c b/src/test_sdis_solve_probe_boundary.c @@ -325,6 +325,8 @@ main(int argc, char** argv) CHK(SOLVE(box_scn, N, iprim, uv, INF, F, 1.0, 0, 0, NULL) == RES_BAD_ARG); CHK(SOLVE(box_scn, N, iprim, uv, INF, F, 1.0, 0, 0, &estimator) == RES_OK); + ref = (H*Tf + LAMBDA * Tb) / (H + LAMBDA); + CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK); CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK); CHK(nfails + nreals == N); @@ -334,12 +336,11 @@ main(int argc, char** argv) CHK(sdis_scene_get_boundary_position(box_scn, iprim, uv, pos) == RES_OK); - ref = (H*Tf + LAMBDA * Tb) / (H + LAMBDA); printf("Boundary 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)); + CHK(eq_eps(T.E, ref, 3*T.SE)); uv[0] = 0.5; iprim = 3; diff --git a/src/test_sdis_volumic_power.c b/src/test_sdis_volumic_power.c @@ -49,6 +49,7 @@ #define T0 320 #define LAMBDA 0.1 #define P0 10 +#define DELTA 1.0/20.0 /******************************************************************************* * Geometry 3D @@ -156,7 +157,7 @@ solid_get_delta { (void)data; CHK(vtx != NULL); - return 1.0/20.0; + return DELTA; } static double diff --git a/src/test_sdis_volumic_power2.c b/src/test_sdis_volumic_power2.c @@ -28,6 +28,16 @@ struct reference { double temperature_3d; /* In celcius */ }; +/* Results in Celcius with delta 0.01 and 10000 realisations + * 0.85: 2D = 190.29 ~ 194.022 +/- 1.87163 [188.407, 199.637]; #failures: 2 + * 0.65: 2D = 259.95 ~ 265.535 +/- 2.14781 [259.091, 271.978]; #failures: 5 + * 0.45: 2D = 286.33 ~ 290.170 +/- 2.25864 [283.394, 296.946]; #failures: 2 + * 0.25: 2D = 235.44 ~ 239.761 +/- 2.27754 [232.928, 246.593]; #failures: 6 + * 0.05: 2D = 192.33 ~ 195.174 +/- 2.22817 [188.489, 201.858]; #failures: 12 + *-0.15; 2D = 156.82 ~ 157.251 +/- 2.13234 [150.854, 163.648]; #failures: 5 + *-0.35; 2D = 123.26 ~ 123.876 +/- 2.00311 [117.866, 129.885]; #failures: 5 + *-0.55; 2D = 90.250 ~ 89.6532 +/- 1.77242 [84.3359, 94.9705]; #failures: 4 */ + static const double vertices[16/*#vertices*/*3/*#coords per vertex*/] = { -0.5,-1.0,-0.5, -0.5, 1.0,-0.5, @@ -234,7 +244,7 @@ main(int argc, char** argv) size_t i; /* In celcius. Computed by EDF with Syrthes */ const struct reference refs[] = { /* Lambda1=1, Lambda2=10, Pw = 10000 */ - {{0, 0.85, 0}, 192.29, 189.13}, + {{0, 0.85, 0}, 190.29, 189.13}, {{0, 0.65, 0}, 259.95, 247.09}, {{0, 0.45, 0}, 286.33, 308.42}, {{0, 0.25, 0}, 235.44, 233.55}, diff --git a/src/test_sdis_volumic_power2_2d.c b/src/test_sdis_volumic_power2_2d.c @@ -25,13 +25,39 @@ #define Tboundary1 NONE #define Tboundary2 NONE #define DELTA 0.01 -#define Tref 286.83 /* In C. Computed with Syrthes at the position 0.5 */ +#define Tref 286.83 /* In Celsius. Computed with Syrthes at the position 0.5 */ /* Dirichlets */ /*#define Tboundary1 373.15*/ /*#define Tboundary2 273.15*/ /*#define DELTA 0.01*/ -/*#define Tref 246.93*/ /* In C. Computed with Syrthes at the position 0.5 */ +/*#define Tref 246.93*/ /* In Celsius. Computed with Syrthes at the position 0.5 */ + +/* >>> Check1: results in Celcius with delta 0.01 and 10000 realisations + * 0.85; 190.29 ~ 192.630 +/- 1.85469 [187.066, 198.194]; #failures: 4 + * 0.65; 259.95 ~ 262.472 +/- 2.20362 [255.861, 269.083]; #failures: 3 + * 0.45; 286.33 ~ 292.721 +/- 2.30152 [285.817, 299.626]; #failures: 3 + * 0.25; 235.44 ~ 238.142 +/- 2.30356 [231.231, 245.052]; #failures: 2 + * 0.05; 192.33 ~ 192.162 +/- 2.18725 [185.601, 198.724]; #failures: 4 + *-0.15; 156.82 ~ 156.178 +/- 2.16574 [149.681, 162.675]; #failures: 8 + *-0.35; 123.26 ~ 123.506 +/- 1.97662 [117.577, 129.436]; #failures: 1 + *-0.55; 90.250 ~ 90.6612 +/- 1.78979 [85.2918, 96.0306]; #failures: 0 + * + * >>> Check2: results in Celcius with delta 0.01 and 10000 realisations + * 0.85; 678.170 ~ 689.735 +/- 13.1487 [650.289, 729.181]; #failures: 12 + * 0.65; 1520.84 ~ 1526.45 +/- 17.2470 [1474.71, 1578.19]; #failures: 38 + * 0.45; 1794.57 ~ 1801.35 +/- 17.5958 [1748.57, 1854.14]; #failures: 31 + * 0.25; 1429.74 ~ 1421.25 +/- 17.4008 [1369.04, 1473.45]; #failures: 32 + * + * >>> Check3: results in Celcius with delta 0.01 and 10000 realisations + * 0.85; 83.99 ~ 84.0036 +/- 0.366646 [82.9037, 85.1035]; #failures: 4 + * 0.65; 73.90 ~ 73.6821 +/- 0.440425 [72.3608, 75.0034]; #failures: 3 + * 0.45; 68.43 ~ 69.4908 +/- 0.460515 [68.1093, 70.8724]; #failures: 3 + * 0.25; 60.61 ~ 60.9922 +/- 0.487816 [59.5287, 62.4556]; #failures: 2 + * 0.05; 52.09 ~ 51.5106 +/- 0.499872 [50.0110, 53.0102]; #failures: 4 + *-0.15; 42.75 ~ 42.1037 +/- 0.493923 [40.6219, 43.5855]; #failures: 8 + *-0.35; 33.04 ~ 33.6234 +/- 0.472444 [32.2060, 35.0407]; #failures: 1 + *-0.55; 24.58 ~ 24.2100 +/- 0.428355 [22.9249, 25.4951]; #failures: 0 */ /* * _\ T1 diff --git a/src/test_sdis_volumic_power4_2d.c b/src/test_sdis_volumic_power4_2d.c @@ -17,13 +17,14 @@ #include "test_sdis_utils.h" #include <rsys/math.h> -#define Tfluid 0 -#define Power 10000 -#define Tboundary -1 -#define Hboundary 50 -#define Lambda 10.0 -#define Delta (1.0/20.0) -#define Nrealisations 10000 +#define Tf1 0 +#define Tf2 100 +#define Power 0 /*10000*/ +#define H1 50 +#define H2 50 +#define LAMBDA 100.0 +#define DELTA (1.0/20.0) +#define N 10000 /* * The 2D scene is a solid slabs stretched along the X dimension to simulate a @@ -47,7 +48,6 @@ * */ - static const double vertices[4/*#vertices*/*2/*#coords per vertex*/] = { -10000.5,-0.5, -10000.5, 0.5, @@ -206,7 +206,8 @@ main(int argc, char** argv) struct interf* interf_param = NULL; struct sdis_device* dev = NULL; struct sdis_data* data = NULL; - struct sdis_medium* fluid = NULL; + struct sdis_medium* fluid1 = NULL; + struct sdis_medium* fluid2 = NULL; struct sdis_medium* solid = NULL; struct sdis_scene* scn = NULL; struct sdis_estimator* estimator = NULL; @@ -214,12 +215,14 @@ main(int argc, char** argv) 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_solid_fluid = NULL; + struct sdis_interface* interf_solid_fluid1 = NULL; + struct sdis_interface* interf_solid_fluid2 = NULL; struct sdis_interface* interfaces[4/*#segment*/]; struct sdis_mc T = SDIS_MC_NULL; + size_t nreals, nfails; double pos[2]; double Tref; - double Tinterf; + double a, b, x; double L; (void)argc, (void)argv; @@ -231,11 +234,19 @@ main(int argc, char** argv) fluid_shader.temperature = fluid_get_temperature; fluid_shader.calorific_capacity = dummy_medium_getter; fluid_shader.volumic_mass = dummy_medium_getter; + + CHK(sdis_data_create + (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data) == RES_OK); + fluid_param = sdis_data_get(data); + fluid_param->temperature = Tf1; + CHK(sdis_fluid_create(dev, &fluid_shader, data, &fluid1) == RES_OK); + CHK(sdis_data_ref_put(data) == RES_OK); + CHK(sdis_data_create (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data) == RES_OK); fluid_param = sdis_data_get(data); - fluid_param->temperature = Tfluid; - CHK(sdis_fluid_create(dev, &fluid_shader, data, &fluid) == RES_OK); + fluid_param->temperature = Tf2; + CHK(sdis_fluid_create(dev, &fluid_shader, data, &fluid2) == RES_OK); CHK(sdis_data_ref_put(data) == RES_OK); /* Setup the solid shader */ @@ -252,8 +263,8 @@ main(int argc, char** argv) solid_param = sdis_data_get(data); solid_param->cp = 500000; solid_param->rho = 1000; - solid_param->lambda = Lambda; - solid_param->delta = Delta; + solid_param->lambda = LAMBDA; + solid_param->delta = DELTA; solid_param->volumic_power = Power; solid_param->temperature = -1; CHK(sdis_solid_create(dev, &solid_shader, data, &solid) == RES_OK); @@ -269,29 +280,40 @@ main(int argc, char** argv) interf_param = sdis_data_get(data); interf_param->h = 0; interf_param->temperature = -1; - CHK(sdis_interface_create(dev, solid, fluid, &interf_shader, data, + CHK(sdis_interface_create(dev, solid, fluid1, &interf_shader, data, &interf_adiabatic) == RES_OK); CHK(sdis_data_ref_put(data) == RES_OK); - /* Create the solid fluid interface */ + /* Create the solid fluid1 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 = H1; + interf_param->temperature = -1; + CHK(sdis_interface_create(dev, solid, fluid1, &interf_shader, data, + &interf_solid_fluid1) == RES_OK); + CHK(sdis_data_ref_put(data) == RES_OK); + + /* Create the solid fluid2 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 = Hboundary; - interf_param->temperature = Tboundary; - CHK(sdis_interface_create(dev, solid, fluid, &interf_shader, data, - &interf_solid_fluid) == RES_OK); + interf_param->h = H2; + interf_param->temperature = -1; + CHK(sdis_interface_create(dev, solid, fluid2, &interf_shader, data, + &interf_solid_fluid2) == 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(fluid1) == RES_OK); + CHK(sdis_medium_ref_put(fluid2) == RES_OK); CHK(sdis_medium_ref_put(solid) == RES_OK); /* Map the interfaces to their square segments */ interfaces[0] = interf_adiabatic; - interfaces[1] = interf_solid_fluid; + interfaces[1] = interf_solid_fluid1; interfaces[2] = interf_adiabatic; - interfaces[3] = interf_solid_fluid; + interfaces[3] = interf_solid_fluid2; #if 0 dump_segments(stdout, vertices, nvertices, indices, nsegments); @@ -304,25 +326,36 @@ main(int argc, char** argv) /* Release the interfaces */ CHK(sdis_interface_ref_put(interf_adiabatic) == RES_OK); - CHK(sdis_interface_ref_put(interf_solid_fluid) == RES_OK); + CHK(sdis_interface_ref_put(interf_solid_fluid1) == RES_OK); + CHK(sdis_interface_ref_put(interf_solid_fluid2) == RES_OK); pos[0] = 0; pos[1] = 0.25; L = vertices[3] - vertices[1]; - if(Tboundary >= 0) { - Tinterf = Tboundary; - } else { - Tinterf = Power*L / (2*Hboundary) + Tfluid; - } - Tref = - Tinterf - + Power / (2*Lambda) * ((L*L)/4.0 - pos[1]*pos[1]); - - CHK(sdis_solve_probe(scn, Nrealisations, pos, INF, 1.f, -1, 0, &estimator) == RES_OK); +#if 1 + x = pos[1] + vertices[3]; + a = (H2*Power*L + H1*H2*(Tf1 - Tf2) + H1*H2*Power*L*L/(2*LAMBDA)) + / (LAMBDA * (H1 + H2) + H1*H2*L); + b = Tf2 + a * LAMBDA / H2; + Tref = -Power / (2*LAMBDA) * x*x + a * x + b; +#else + tmp = LAMBDA / L; + T1 = H1 * (H2+tmp) / (tmp*(H1+H2) + H1*H2) * Tf1 + + H2 * tmp / (tmp*(H1+H2) + H1*H2) * Tf2; + T2 = H1 * tmp / (tmp*(H1+H2) + H1*H2) * Tf1 + + H2 * (H1+tmp) / (tmp*(H1+H2) + H1*H2) * Tf2; + Tref = T2 + (T1-T2)/L * (pos[1] + vertices[3]); +#endif + + 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_get_realisation_count(estimator, &nreals) == RES_OK); + CHK(sdis_estimator_get_failure_count(estimator, &nfails) == 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); + printf("#realisations: %lu; #failures: %lu\n", + (unsigned long)nreals, (unsigned long)nfails); CHK(eq_eps(T.E, Tref, T.SE*3)); CHK(sdis_estimator_ref_put(estimator) == RES_OK);