stardis-solver

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

commit 1c05f25347f2c49cb463a6025921e8b63bb124a5
parent 49760dd523b8cfb1b779a047dd26c5b7ca0899fb
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date:   Tue,  1 Sep 2020 18:20:19 +0200

Add tests on green

Diffstat:
Msrc/test_sdis_conducto_radiative.c | 49+++++++++++++++++++++++++++++++++++++++++++------
Msrc/test_sdis_conducto_radiative_2d.c | 2+-
Msrc/test_sdis_flux.c | 83++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---------------
Msrc/test_sdis_solve_probe.c | 30++++++++++++++++++++++++++++--
Msrc/test_sdis_volumic_power.c | 102++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---------------------
5 files changed, 215 insertions(+), 51 deletions(-)

diff --git a/src/test_sdis_conducto_radiative.c b/src/test_sdis_conducto_radiative.c @@ -278,6 +278,7 @@ main(int argc, char** argv) const double T1 = 310; /* Fixed temperature on the right side of the system */ const double thickness = 2.0; /* Thickness of the solid along X */ double Ts0, Ts1, hr, tmp; + struct interfac* p_intface; (void)argc, (void)argv; OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator)); @@ -376,11 +377,10 @@ main(int argc, char** argv) get_position, &geom, &scn)); hr = 4.0 * BOLTZMANN_CONSTANT * Tref*Tref*Tref * emissivity; - tmp = lambda/(2*lambda + thickness*hr) * (T1 - T0); - Ts0 = T0 + tmp; - Ts1 = T1 - tmp; /* Run the simulations */ + p_intface + = (struct interfac*)sdis_data_get(sdis_interface_get_data(interfaces[4])); OK(ssp_rng_create(&allocator, &ssp_rng_kiss, &rng)); FOR_EACH(isimul, 0, nsimuls) { struct sdis_mc T = SDIS_MC_NULL; @@ -393,6 +393,10 @@ main(int argc, char** argv) size_t nreals = 0; size_t nfails = 0; const size_t N = 10000; + double T1b; + + /* Reset temperature */ + p_intface->temperature = T1; solve_args.nrealisations = N; solve_args.time_range[0] = INF; @@ -400,6 +404,7 @@ main(int argc, char** argv) solve_args.position[0] = ssp_rng_uniform_double(rng, -0.9, 0.9); solve_args.position[1] = ssp_rng_uniform_double(rng, -0.9, 0.9); solve_args.position[2] = ssp_rng_uniform_double(rng, -0.9, 0.9); + u = (solve_args.position[0] + 1) / thickness; solve_args.reference_temperature = Tref; OK(sdis_solve_probe(scn, &solve_args, &estimator)); @@ -408,10 +413,12 @@ main(int argc, char** argv) OK(sdis_estimator_get_temperature(estimator, &T)); OK(sdis_estimator_get_realisation_time(estimator, &time)); - u = (solve_args.position[0] + 1) / thickness; + tmp = lambda / (2 * lambda + thickness * hr) * (T1 - T0); + Ts0 = T0 + tmp; + Ts1 = T1 - tmp; ref = u * Ts1 + (1-u) * Ts0; - printf("Temperature at (%g, %g, %g) = %g ~ %g +/- %g\n", - SPLIT3(solve_args.position), ref, T.E, T.SE); + printf("Temperature at (%g, %g, %g) with T1=%g = %g ~ %g +/- %g\n", + SPLIT3(solve_args.position), p_intface->temperature, ref, T.E, T.SE); printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); @@ -427,6 +434,36 @@ main(int argc, char** argv) OK(sdis_estimator_ref_put(estimator)); OK(sdis_estimator_ref_put(estimator2)); + printf("\n"); + + /* Check green used at a different temperature */ + p_intface->temperature = T1b = T1 + (isimul + 1) * 10; + + OK(sdis_solve_probe(scn, &solve_args, &estimator)); + OK(sdis_estimator_get_realisation_count(estimator, &nreals)); + OK(sdis_estimator_get_failure_count(estimator, &nfails)); + OK(sdis_estimator_get_temperature(estimator, &T)); + OK(sdis_estimator_get_realisation_time(estimator, &time)); + + tmp = lambda / (2 * lambda + thickness * hr) * (T1b - T0); + Ts0 = T0 + tmp; + Ts1 = T1b - tmp; + ref = u * Ts1 + (1 - u) * Ts0; + printf("Temperature at (%g, %g, %g) with T1=%g = %g ~ %g +/- %g\n", + SPLIT3(solve_args.position), p_intface->temperature, ref, T.E, T.SE); + printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); + printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); + + CHK(nfails + nreals == N); + CHK(nfails < N / 1000); + CHK(eq_eps(T.E, ref, 3*T.SE) == 1); + + OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + check_green_function(green); + check_estimator_eq(estimator, estimator2); + + OK(sdis_estimator_ref_put(estimator)); + OK(sdis_estimator_ref_put(estimator2)); OK(sdis_green_function_ref_put(green)); solve_args.nrealisations = 10; diff --git a/src/test_sdis_conducto_radiative_2d.c b/src/test_sdis_conducto_radiative_2d.c @@ -415,7 +415,7 @@ main(int argc, char** argv) u = (solve_args.position[0] + 1) / thickness; ref = u * Ts1 + (1-u) * Ts0; - printf("Temperature at (%g, %g) = %g ~ %g +/- %g\n", + printf("Temperature at (%g, %g) = %g ~ %g +/- %g\n", SPLIT2(solve_args.position), ref, T.E, T.SE); printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); diff --git a/src/test_sdis_flux.c b/src/test_sdis_flux.c @@ -127,7 +127,10 @@ interface_get_flux * Helper functions ******************************************************************************/ static void -solve(struct sdis_scene* scn, const double pos[]) +solve + (struct sdis_scene* scn, + const double pos[], + struct interf* interf) { char dump[128]; struct time t0, t1, t2; @@ -142,9 +145,12 @@ solve(struct sdis_scene* scn, const double pos[]) double ref; const double time_range[2] = {INF, INF}; enum sdis_scene_dimension dim; - ASSERT(scn && pos); + ASSERT(scn && pos && interf); - ref = T0 + (1 - pos[0]) * PHI/LAMBDA; + /* Restore phi value */ + interf->phi = PHI; + + ref = T0 + (1 - pos[0]) * interf->phi/LAMBDA; OK(sdis_scene_get_dimension(scn, &dim)); @@ -167,18 +173,18 @@ solve(struct sdis_scene* scn, const double pos[]) switch(dim) { case SDIS_SCENE_2D: - printf("Temperature at (%g %g) = %g ~ %g +/- %g\n", - SPLIT2(pos), ref, T.E, T.SE); + printf("Temperature at (%g %g) with phi=%g = %g ~ %g +/- %g\n", + SPLIT2(pos), interf->phi, ref, T.E, T.SE); break; case SDIS_SCENE_3D: - printf("Temperature at (%g %g %g) = %g ~ %g +/- %g\n", - SPLIT3(pos), ref, T.E, T.SE); + printf("Temperature at (%g %g %g) with phi=%g = %g ~ %g +/- %g\n", + SPLIT3(pos), interf->phi, ref, T.E, T.SE); break; default: FATAL("Unreachable code.\n"); break; } - printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); - printf("Elapsed time = %s\n\n", dump); + printf("Elapsed time = %s\n", dump); + printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE); CHK(nfails + nreals == N); CHK(nfails < N/1000); @@ -196,12 +202,12 @@ solve(struct sdis_scene* scn, const double pos[]) switch(dim) { case SDIS_SCENE_2D: - printf("Green temperature at (%g %g) = %g ~ %g +/- %g\n", - SPLIT2(pos), ref, T.E, T.SE); + printf("Green temperature at (%g %g) with phi=%g = %g ~ %g +/- %g\n", + SPLIT2(pos), interf->phi, ref, T.E, T.SE); break; case SDIS_SCENE_3D: - printf("Green temperature at (%g %g %g) = %g ~ %g +/- %g\n", - SPLIT3(pos), ref, T.E, T.SE); + printf("Green temperature at (%g %g %g) with phi=%g = %g ~ %g +/- %g\n", + SPLIT3(pos), interf->phi, ref, T.E, T.SE); break; default: FATAL("Unreachable code.\n"); break; } @@ -218,6 +224,53 @@ solve(struct sdis_scene* scn, const double pos[]) OK(sdis_estimator_ref_put(estimator)); OK(sdis_estimator_ref_put(estimator2)); + printf("\n"); + + /* Check green used at a different phi */ + interf->phi = 3 * PHI; + + time_current(&t0); + OK(sdis_solve_probe(scn, &solve_args, &estimator)); + time_sub(&t0, time_current(&t1), &t0); + time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); + + OK(sdis_estimator_get_realisation_count(estimator, &nreals)); + OK(sdis_estimator_get_failure_count(estimator, &nfails)); + OK(sdis_estimator_get_temperature(estimator, &T)); + OK(sdis_estimator_get_realisation_time(estimator, &time)); + + ref = T0 + (1 - pos[0]) * interf->phi/LAMBDA; + + switch (dim) { + case SDIS_SCENE_2D: + printf("Temperature at (%g %g) with phi=%g = %g ~ %g +/- %g\n", + SPLIT2(pos), interf->phi, ref, T.E, T.SE); + break; + case SDIS_SCENE_3D: + printf("Temperature at (%g %g %g) with phi=%g = %g ~ %g +/- %g\n", + SPLIT3(pos), interf->phi, ref, T.E, T.SE); + break; + default: FATAL("Unreachable code.\n"); break; + } + printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); + printf("Elapsed time = %s\n", dump); + printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE); + + CHK(nfails + nreals == N); + CHK(nfails < N / 1000); + CHK(eq_eps(T.E, ref, T.SE * 3)); + + time_current(&t0); + OK(sdis_green_function_solve(green, time_range, &estimator2)); + time_sub(&t0, time_current(&t1), &t0); + time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); + printf("Green solve time = %s\n", dump); + + check_green_function(green); + check_estimator_eq(estimator, estimator2); + + OK(sdis_estimator_ref_put(estimator)); + OK(sdis_estimator_ref_put(estimator2)); OK(sdis_green_function_ref_put(green)); printf("\n"); @@ -328,9 +381,9 @@ main(int argc, char** argv) /* Solve */ d3_splat(pos, 0.25); printf(">> Box scene\n"); - solve(box_scn, pos); + solve(box_scn, pos, interf_props); printf(">> Square Scene\n"); - solve(square_scn, pos); + solve(square_scn, pos, interf_props); OK(sdis_scene_ref_put(box_scn)); OK(sdis_scene_ref_put(square_scn)); diff --git a/src/test_sdis_solve_probe.c b/src/test_sdis_solve_probe.c @@ -403,8 +403,8 @@ main(int argc, char** argv) OK(sdis_estimator_get_rng_state(estimator, &rng_state)); ref = 300; - printf("Temperature at (%g, %g, %g) = %g ~ %g +/- %g\n", - SPLIT3(solve_args.position), ref, T.E, T.SE); + printf("Temperature at (%g, %g, %g) with Tfluid=%g = %g ~ %g +/- %g\n", + SPLIT3(solve_args.position), fluid_param->temperature, ref, T.E, T.SE); printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); @@ -444,6 +444,32 @@ main(int argc, char** argv) check_green_function(green); check_estimator_eq(estimator, estimator2); + OK(sdis_estimator_ref_put(estimator)); + OK(sdis_estimator_ref_put(estimator2)); + printf("\n"); + + /* Check green used at a different temperature */ + fluid_param->temperature = 500; + OK(sdis_solve_probe(scn, &solve_args, &estimator)); + OK(sdis_estimator_get_realisation_count(estimator, &nreals)); + OK(sdis_estimator_get_failure_count(estimator, &nfails)); + OK(sdis_estimator_get_temperature(estimator, &T)); + OK(sdis_estimator_get_realisation_time(estimator, &time)); + + ref = 500; + printf("Temperature at (%g, %g, %g) with Tfluid=%g = %g ~ %g +/- %g\n", + SPLIT3(solve_args.position), fluid_param->temperature, ref, T.E, T.SE); + printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); + printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); + + CHK(nfails + nreals == N); + CHK(nfails < N / 1000); + CHK(eq_eps(T.E, ref, T.SE)); + + OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + check_green_function(green); + check_estimator_eq(estimator, estimator2); + BA(sdis_green_function_ref_get(NULL)); OK(sdis_green_function_ref_get(green)); BA(sdis_green_function_ref_put(NULL)); diff --git a/src/test_sdis_volumic_power.c b/src/test_sdis_volumic_power.c @@ -60,6 +60,7 @@ struct solid { double rho; double cp; double delta; + double vpower; }; static double @@ -118,7 +119,7 @@ solid_get_volumic_power { (void)data; CHK(vtx != NULL); - return P0; + return ((struct solid*)sdis_data_cget(data))->vpower; } /******************************************************************************* @@ -149,7 +150,10 @@ interface_get_convection_coef * Helper functions ******************************************************************************/ static void -solve(struct sdis_scene* scn, const double pos[]) +solve + (struct sdis_scene* scn, + const double pos[], + struct solid* solid) { char dump[128]; struct time t0, t1, t2; @@ -157,17 +161,21 @@ solve(struct sdis_scene* scn, const double pos[]) struct sdis_estimator* estimator2; struct sdis_green_function* green; struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; - struct sdis_mc T; + struct sdis_mc T = SDIS_MC_NULL; + struct sdis_mc time = SDIS_MC_NULL; size_t nreals; size_t nfails; double x; double ref; const double time_range[2] = {INF, INF}; enum sdis_scene_dimension dim; - ASSERT(scn && pos); + ASSERT(scn && pos && solid); + + /* Restore power value */ + solid->vpower = P0; x = pos[0] - 0.5; - ref = P0 / (2*LAMBDA) * (1.0/4.0 - x*x) + T0; + ref = solid->vpower / (2*LAMBDA) * (1.0/4.0 - x*x) + T0; OK(sdis_scene_get_dimension(scn, &dim)); @@ -186,25 +194,28 @@ solve(struct sdis_scene* scn, const double pos[]) OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); OK(sdis_estimator_get_temperature(estimator, &T)); + OK(sdis_estimator_get_realisation_time(estimator, &time)); switch(dim) { case SDIS_SCENE_2D: - printf("Temperature at (%g %g) = %g ~ %g +/- %g [%g, %g]\n", - SPLIT2(pos), ref, T.E, T.SE, T.E-3*T.SE, T.E+3*T.SE); + printf("Temperature at (%g %g) with Power=%g = %g ~ %g +/- %g [%g, %g]\n", + SPLIT2(pos), solid->vpower, ref, T.E, T.SE, T.E-3*T.SE, T.E+3*T.SE); break; case SDIS_SCENE_3D: - printf("Temperature at (%g %g %g) = %g ~ %g +/- %g [%g, %g]\n", - SPLIT3(pos), ref, T.E, T.SE, T.E -3*T.SE, T.E + 3*T.SE); + printf("Temperature at (%g %g %g)th Power=%g = %g ~ %g +/- %g [%g, %g]\n", + SPLIT3(pos), solid->vpower, ref, T.E, T.SE, T.E -3*T.SE, T.E + 3*T.SE); break; default: FATAL("Unreachable code.\n"); break; } printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); - printf("Elapsed time = %s\n\n", dump); + printf("Elapsed time = %s\n", dump); + printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE); CHK(nfails + nreals == N); CHK(nfails < N/1000); CHK(eq_eps(T.E, ref, T.SE*3)); + /* Check green function */ time_current(&t0); OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); time_current(&t1); @@ -217,12 +228,12 @@ solve(struct sdis_scene* scn, const double pos[]) switch(dim) { case SDIS_SCENE_2D: - printf("Green temperature at (%g %g) = %g ~ %g +/- %g\n", - SPLIT2(pos), ref, T.E, T.SE); + printf("Green temperature at (%g %g) with Power=%g = %g ~ %g +/- %g\n", + SPLIT2(pos), solid->vpower, ref, T.E, T.SE); break; case SDIS_SCENE_3D: - printf("Green temperature at (%g %g %g) = %g ~ %g +/- %g\n", - SPLIT3(pos), ref, T.E, T.SE); + printf("Green temperature at (%g %g %g) with Power=%g = %g ~ %g +/- %g\n", + SPLIT3(pos), solid->vpower, ref, T.E, T.SE); break; default: FATAL("Unreachable code.\n"); break; } @@ -239,6 +250,53 @@ solve(struct sdis_scene* scn, const double pos[]) OK(sdis_estimator_ref_put(estimator)); OK(sdis_estimator_ref_put(estimator2)); + printf("\n"); + + /* Check green used at a different power */ + solid->vpower = 3 * P0; + + time_current(&t0); + OK(sdis_solve_probe(scn, &solve_args, &estimator)); + time_sub(&t0, time_current(&t1), &t0); + time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); + + OK(sdis_estimator_get_realisation_count(estimator, &nreals)); + OK(sdis_estimator_get_failure_count(estimator, &nfails)); + OK(sdis_estimator_get_temperature(estimator, &T)); + OK(sdis_estimator_get_realisation_time(estimator, &time)); + + ref = solid->vpower / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0; + + switch (dim) { + case SDIS_SCENE_2D: + printf("Temperature at (%g %g) with Power=%g = %g ~ %g +/- %g [%g, %g]\n", + SPLIT2(pos), solid->vpower, ref, T.E, T.SE, T.E - 3 * T.SE, T.E + 3 * T.SE); + break; + case SDIS_SCENE_3D: + printf("Temperature at (%g %g %g) with Power=%g = %g ~ %g +/- %g [%g, %g]\n", + SPLIT3(pos), solid->vpower, ref, T.E, T.SE, T.E - 3 * T.SE, T.E + 3 * T.SE); + break; + default: FATAL("Unreachable code.\n"); break; + } + printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); + printf("Elapsed time = %s\n", dump); + printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); + + CHK(nfails + nreals == N); + CHK(nfails < N / 1000); + CHK(eq_eps(T.E, ref, T.SE * 3)); + + time_current(&t0); + OK(sdis_green_function_solve(green, time_range, &estimator2)); + time_sub(&t0, time_current(&t1), &t0); + time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); + printf("Green solve time = %s\n", dump); + + check_green_function(green); + check_estimator_eq(estimator, estimator2); + + OK(sdis_estimator_ref_put(estimator)); + OK(sdis_estimator_ref_put(estimator2)); OK(sdis_green_function_ref_put(green)); printf("\n"); @@ -255,7 +313,6 @@ main(int argc, char** argv) struct sdis_device* dev = NULL; struct sdis_medium* fluid = NULL; struct sdis_medium* solid = NULL; - struct sdis_medium* solid2 = NULL; /* For debug */ struct sdis_interface* interf_adiabatic = NULL; struct sdis_interface* interf_T0 = NULL; struct sdis_scene* box_scn = NULL; @@ -291,18 +348,10 @@ main(int argc, char** argv) solid_props->cp = 2; solid_props->rho = 25; solid_props->delta = DELTA; + solid_props->vpower = P0; OK(sdis_solid_create(dev, &solid_shader, data, &solid)); OK(sdis_data_ref_put(data)); - OK(sdis_data_create(dev, sizeof(struct solid), 16, NULL, &data)); - solid_props = sdis_data_get(data); - solid_props->lambda = 0; - solid_props->cp = 0; - solid_props->rho = 0; - solid_props->delta = DELTA/4; - OK(sdis_solid_create(dev, &solid_shader, data, &solid2)); - OK(sdis_data_ref_put(data)); - /* Setup the interface shader */ interf_shader.convection_coef = interface_get_convection_coef; interf_shader.front.temperature = interface_get_temperature; @@ -325,7 +374,6 @@ main(int argc, char** argv) /* Release the media */ OK(sdis_medium_ref_put(solid)); - OK(sdis_medium_ref_put(solid2)); OK(sdis_medium_ref_put(fluid)); /* Map the interfaces to their box triangles */ @@ -359,9 +407,9 @@ main(int argc, char** argv) /* Solve */ d3_splat(pos, 0.25); printf(">> Box scene\n"); - solve(box_scn, pos); + solve(box_scn, pos, solid_props); printf(">> Square scene\n"); - solve(square_scn, pos); + solve(square_scn, pos, solid_props); OK(sdis_scene_ref_put(box_scn)); OK(sdis_scene_ref_put(square_scn));