commit a5af128204148e01af2975d223aeec8771aaa71c
parent 68995e90d8229cc41bac349df4ecfb33d375c221
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Tue, 8 Dec 2020 11:03:48 +0100
Fix flux computation on surfaces (was computing a density)
Diffstat:
5 files changed, 94 insertions(+), 48 deletions(-)
diff --git a/doc/stardis-output.5.txt b/doc/stardis-output.5.txt
@@ -139,34 +139,36 @@ _______
<failures-report> ::= <error-count> <success-count>
<ext-probe-temp> ::= "Temperature at" <ext-probe-position> <ext-time> \
- <ext-MC-estimate>
+ <ext-temp-estimate>
<ext-failures-report>
<ext-mean-temp> ::= "Temperature at boundary" <file-name> <ext-time> \
- <ext-MC-estimate>
+ <ext-temp-estimate>
<ext-failures-report>
<ext-mean-flux> ::= "Temperature at boundary" <file-name> <ext-time> \
- <ext-MC-estimate>
+ <ext-temp-estimate>
"Convective flux at boundary " <file-name> \
- <ext-time> <ext-MC-estimate>
+ <ext-time> <ext-flux-estimate>
"Radiative flux at boundary" <file-name> <ext-time> \
- <ext-MC-estimate>
+ <ext-flux-estimate>
"Imposed flux at boundary" <file-name> <ext-time> \
- <ext-MC-estimate>
+ <ext-flux-estimate>
"Total flux at boundary" <file-name> <ext-time> \
- <ext-MC-estimate>
+ <ext-flux-estimate>
<ext-failures-report>
<ext-medium-temp> ::= "Temperature in medium" <medium-name> <ext-time> \
- <ext-MC-estimate>
+ <ext-temp-estimate>
<ext-failures-report>
<ext-probe-position> ::= "[" <x-value> "," <y-value> "," <z-value> "]"
<ext-time> ::= "at t=" <time-value> "="
-<ext-MC-estimate> ::= <expected-value> "+/-" <standard-deviation>
+<ext-temp-estimate> ::= <expected-value> "K +/-" <standard-deviation>
+
+<ext-flux-estimate> ::= <expected-value> "W +/-" <standard-deviation>
<ext-failures-report> ::= "#failures:" <error-count> "/" <success-count>
diff --git a/src/stardis-app.c b/src/stardis-app.c
@@ -229,6 +229,7 @@ stardis_init
darray_size_t_init(stardis->allocator, &stardis->compute_surface.primitives);
darray_sides_init(stardis->allocator, &stardis->compute_surface.sides);
darray_uint_init(stardis->allocator, &stardis->compute_surface.err_triangles);
+ stardis->compute_surface.area = 0;
stardis->samples = args->samples;
stardis->nthreads = args->nthreads;
stardis->scale_factor = -1; /* invalid value */
@@ -309,6 +310,8 @@ stardis_init
goto error;
}
}
+ logger_print(stardis->logger, LOG_OUTPUT, "Compute surface area is %g m2\n",
+ stardis->compute_surface.area);
}
/* Create enclosures */
diff --git a/src/stardis-app.h b/src/stardis-app.h
@@ -712,6 +712,7 @@ struct compute_surface {
struct darray_size_t primitives;
struct darray_sides sides;
struct darray_uint err_triangles;
+ double area; /* in m2, regardless of scale factor */
};
struct stardis {
diff --git a/src/stardis-compute.c b/src/stardis-compute.c
@@ -972,6 +972,8 @@ read_compute_surface
struct add_geom_ctx add_geom_ctx;
struct sg3d_geometry_add_callbacks callbacks = SG3D_ADD_CALLBACKS_NULL__;
const char* file;
+ size_t i;
+ double _2area = 0;
ASSERT(stardis);
@@ -1006,6 +1008,26 @@ read_compute_surface
goto error;
}
+ /* Compute area of the compute surface */
+ FOR_EACH(i, 0, add_geom_ctx.stl_desc.triangles_count) {
+ const unsigned* trg = add_geom_ctx.stl_desc.indices + (i * 3);
+ const float* v0f = add_geom_ctx.stl_desc.vertices + (trg[0] * 3);
+ const float* v1f = add_geom_ctx.stl_desc.vertices + (trg[1] * 3);
+ const float* v2f = add_geom_ctx.stl_desc.vertices + (trg[2] * 3);
+ double v0[3], v1[3], v2[3], edge0[3], edge1[3], normal[3], norm;
+
+ /* Compute component area and volume */
+ d3_sub(edge0, d3_set_f3(v1, v1f), d3_set_f3(v0, v0f));
+ d3_sub(edge1, d3_set_f3(v2, v2f), d3_set_f3(v0, v0f));
+ d3_cross(normal, edge0, edge1);
+ norm = d3_normalize(normal, normal);
+ ASSERT(norm);
+ _2area += norm;
+ }
+
+ stardis->compute_surface.area = _2area * 0.5
+ * stardis->scale_factor * stardis->scale_factor;
+
end:
if(sstl) SSTL(ref_put(sstl));
return res;
diff --git a/src/stardis-output.c b/src/stardis-output.c
@@ -1330,12 +1330,13 @@ print_single_MC_result
case MODE_PROBE_COMPUTE:
if(stardis->mode & MODE_EXTENDED_RESULTS) {
if(stardis->time_range[0] == stardis->time_range[1])
- fprintf(stream, "Temperature at [%g, %g, %g] at t=%g = %g +/- %g\n",
+ fprintf(stream, "Temperature at [%g, %g, %g] at t=%g = %g K +/- %g\n",
SPLIT3(stardis->probe), stardis->time_range[0],
result.E, /* Expected value */
result.SE); /* Standard error */
else
- fprintf(stream, "Temperature at [%g, %g, %g] with t in [%g %g] = %g +/- %g\n",
+ fprintf(stream,
+ "Temperature at [%g, %g, %g] with t in [%g %g] = %g K +/- %g\n",
SPLIT3(stardis->probe), SPLIT2(stardis->time_range),
result.E, /* Expected value */
result.SE); /* Standard error */
@@ -1346,12 +1347,14 @@ print_single_MC_result
case MODE_PROBE_COMPUTE_ON_INTERFACE:
if(stardis->mode & MODE_EXTENDED_RESULTS) {
if(stardis->time_range[0] == stardis->time_range[1])
- fprintf(stream, "Boundary temperature at [%g, %g, %g] at t=%g = %g +/- %g\n",
+ fprintf(stream,
+ "Boundary temperature at [%g, %g, %g] at t=%g = %g K +/- %g\n",
SPLIT3(stardis->probe), stardis->time_range[0],
result.E, /* Expected value */
result.SE); /* Standard error */
else
- fprintf(stream, "Boundary temperature at [%g, %g, %g] with t in [%g %g] = %g +/- %g\n",
+ fprintf(stream,
+ "Boundary temperature at [%g, %g, %g] with t in [%g %g] = %g K +/- %g\n",
SPLIT3(stardis->probe), SPLIT2(stardis->time_range),
result.E, /* Expected value */
result.SE); /* Standard error */
@@ -1362,12 +1365,13 @@ print_single_MC_result
case MODE_MEDIUM_COMPUTE:
if(stardis->mode & MODE_EXTENDED_RESULTS) {
if(stardis->time_range[0] == stardis->time_range[1])
- fprintf(stream, "Temperature in medium '%s' at t=%g = %g +/- %g\n",
+ fprintf(stream, "Temperature in medium '%s' at t=%g = %g K +/- %g\n",
str_cget(&stardis->solve_name), stardis->time_range[0],
result.E, /* Expected value */
result.SE); /* Standard error */
else
- fprintf(stream, "Temperature in medium '%s' with t in [%g %g] = %g +/- %g\n",
+ fprintf(stream,
+ "Temperature in medium '%s' with t in [%g %g] = %g K +/- %g\n",
str_cget(&stardis->solve_name), SPLIT2(stardis->time_range),
result.E, /* Expected value */
result.SE); /* Standard error */
@@ -1378,12 +1382,13 @@ print_single_MC_result
case MODE_BOUNDARY_COMPUTE:
if(stardis->mode & MODE_EXTENDED_RESULTS) {
if(stardis->time_range[0] == stardis->time_range[1])
- fprintf(stream, "Temperature at boundary '%s' at t=%g = %g +/- %g\n",
+ fprintf(stream, "Temperature at boundary '%s' at t=%g = %g K +/- %g\n",
str_cget(&stardis->solve_name), stardis->time_range[0],
result.E, /* Expected value */
result.SE); /* Standard error */
else
- fprintf(stream, "Temperature at boundary '%s' with t in [%g %g] = %g +/- %g\n",
+ fprintf(stream,
+ "Temperature at boundary '%s' with t in [%g %g] = %g K +/- %g\n",
str_cget(&stardis->solve_name), SPLIT2(stardis->time_range),
result.E, /* Expected value */
result.SE); /* Standard error */
@@ -1398,69 +1403,82 @@ print_single_MC_result
if(stardis->mode & MODE_EXTENDED_RESULTS) {
if(stardis->time_range[0] == stardis->time_range[1])
- fprintf(stream, "Temperature at boundary '%s' at t=%g = %g +/- %g\n",
+ fprintf(stream, "Temperature at boundary '%s' at t=%g = %g K +/- %g\n",
str_cget(&stardis->solve_name), stardis->time_range[0],
result.E, /* Expected value */
result.SE); /* Standard error */
else
- fprintf(stream, "Temperature at boundary '%s' with t in [%g %g] = %g +/- %g\n",
+ fprintf(stream,
+ "Temperature at boundary '%s' with t in [%g %g] = %g K +/- %g\n",
str_cget(&stardis->solve_name), SPLIT2(stardis->time_range),
result.E, /* Expected value */
result.SE); /* Standard error */
ERR(sdis_estimator_get_convective_flux(estimator, &result));
if(stardis->time_range[0] == stardis->time_range[1])
- fprintf(stream, "Convective flux at boundary '%s' at t=%g = %g +/- %g\n",
+ fprintf(stream, "Convective flux at boundary '%s' at t=%g = %g W +/- %g\n",
str_cget(&stardis->solve_name), stardis->time_range[0],
- result.E, /* Expected value */
- result.SE); /* Standard error */
+ stardis->compute_surface.area * result.E, /* Expected value */
+ stardis->compute_surface.area * result.SE); /* Standard error */
else
- fprintf(stream, "Convective flux at boundary '%s' with t in [%g %g] = %g +/- %g\n",
+ fprintf(stream,
+ "Convective flux at boundary '%s' with t in [%g %g] = %g W +/- %g\n",
str_cget(&stardis->solve_name), SPLIT2(stardis->time_range),
- result.E, /* Expected value */
- result.SE); /* Standard error */
+ stardis->compute_surface.area * result.E, /* Expected value */
+ stardis->compute_surface.area * result.SE); /* Standard error */
ERR(sdis_estimator_get_radiative_flux(estimator, &result));
if(stardis->time_range[0] == stardis->time_range[1])
- fprintf(stream, "Radiative flux at boundary '%s' at t=%g = %g +/- %g\n",
+ fprintf(stream, "Radiative flux at boundary '%s' at t=%g = %g W +/- %g\n",
str_cget(&stardis->solve_name), stardis->time_range[0],
- result.E, /* Expected value */
- result.SE); /* Standard error */
+ stardis->compute_surface.area * result.E, /* Expected value */
+ stardis->compute_surface.area * result.SE); /* Standard error */
else
- fprintf(stream, "Radiative flux at boundary '%s' with t in [%g %g] = %g +/- %g\n",
+ fprintf(stream,
+ "Radiative flux at boundary '%s' with t in [%g %g] = %g W +/- %g\n",
str_cget(&stardis->solve_name), SPLIT2(stardis->time_range),
- result.E, /* Expected value */
- result.SE); /* Standard error */
+ stardis->compute_surface.area * result.E, /* Expected value */
+ stardis->compute_surface.area * result.SE); /* Standard error */
ERR(sdis_estimator_get_imposed_flux(estimator, &result));
if(stardis->time_range[0] == stardis->time_range[1])
- fprintf(stream, "Imposed flux at boundary '%s' at t=%g = %g +/- %g\n",
+ fprintf(stream, "Imposed flux at boundary '%s' at t=%g = %g W +/- %g\n",
str_cget(&stardis->solve_name), stardis->time_range[0],
- result.E, /* Expected value */
- result.SE); /* Standard error */
+ stardis->compute_surface.area * result.E, /* Expected value */
+ stardis->compute_surface.area * result.SE); /* Standard error */
else
- fprintf(stream, "Imposed flux at boundary '%s' with t in [%g %g] = %g +/- %g\n",
+ fprintf(stream,
+ "Imposed flux at boundary '%s' with t in [%g %g] = %g W +/- %g\n",
str_cget(&stardis->solve_name), SPLIT2(stardis->time_range),
- result.E, /* Expected value */
- result.SE); /* Standard error */
+ stardis->compute_surface.area * result.E, /* Expected value */
+ stardis->compute_surface.area * result.SE); /* Standard error */
ERR(sdis_estimator_get_total_flux(estimator, &result));
if(stardis->time_range[0] == stardis->time_range[1])
- fprintf(stream, "Total flux Flux at boundary '%s' at t=%g = %g +/- %g\n",
+ fprintf(stream, "Total flux Flux at boundary '%s' at t=%g = %g W +/- %g\n",
str_cget(&stardis->solve_name), stardis->time_range[0],
- result.E, /* Expected value */
- result.SE); /* Standard error */
+ stardis->compute_surface.area * result.E, /* Expected value */
+ stardis->compute_surface.area * result.SE); /* Standard error */
else
- fprintf(stream, "Total flux Flux at boundary '%s' with t in [%g %g] = %g +/- %g\n",
+ fprintf(stream,
+ "Total flux Flux at boundary '%s' with t in [%g %g] = %g W +/- %g\n",
str_cget(&stardis->solve_name), SPLIT2(stardis->time_range),
- result.E, /* Expected value */
- result.SE); /* Standard error */
+ stardis->compute_surface.area * result.E, /* Expected value */
+ stardis->compute_surface.area * result.SE); /* Standard error */
} else {
- fprintf(stream, "%g %g ", result.E, result.SE);
+ fprintf(stream, "%g %g ", result.E, result.SE); /* Temperature */
ERR(sdis_estimator_get_convective_flux(estimator, &result));
- fprintf(stream, "%g %g ", result.E, result.SE);
+ fprintf(stream, "%g %g ",
+ stardis->compute_surface.area * result.E,
+ stardis->compute_surface.area * result.SE);
ERR(sdis_estimator_get_radiative_flux(estimator, &result));
- fprintf(stream, "%g %g ", result.E, result.SE);
+ fprintf(stream, "%g %g ",
+ stardis->compute_surface.area * result.E,
+ stardis->compute_surface.area * result.SE);
ERR(sdis_estimator_get_imposed_flux(estimator, &result));
- fprintf(stream, "%g %g ", result.E, result.SE);
+ fprintf(stream, "%g %g ",
+ stardis->compute_surface.area * result.E,
+ stardis->compute_surface.area * result.SE);
ERR(sdis_estimator_get_total_flux(estimator, &result));
- fprintf(stream, "%g %g ", result.E, result.SE);
+ fprintf(stream, "%g %g ",
+ stardis->compute_surface.area * result.E,
+ stardis->compute_surface.area * result.SE);
fprintf(stream, "%lu %lu\n", nfailures, nsamples);
}
break;
