commit c6b6ba86fe0f596505abc6787e9ba7371d409c37
parent 63bfef8cfe6a92bc2605cea3ed5689ac67f78262
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Tue, 29 Sep 2020 14:08:54 +0200
Add imposed_flux management to solve_boundary_flux API calls
Improve detection of invalid cases too
Diffstat:
5 files changed, 100 insertions(+), 10 deletions(-)
diff --git a/src/sdis.h b/src/sdis.h
@@ -940,6 +940,11 @@ sdis_estimator_get_radiative_flux
struct sdis_mc* flux);
SDIS_API res_T
+sdis_estimator_get_imposed_flux
+ (const struct sdis_estimator* estimator,
+ struct sdis_mc* flux);
+
+SDIS_API res_T
sdis_estimator_get_total_flux
(const struct sdis_estimator* estimator,
struct sdis_mc* flux);
diff --git a/src/sdis_estimator.c b/src/sdis_estimator.c
@@ -136,6 +136,16 @@ sdis_estimator_get_radiative_flux
}
res_T
+sdis_estimator_get_imposed_flux
+ (const struct sdis_estimator* estimator, struct sdis_mc* flux)
+{
+ if(!estimator || !flux || estimator->type != SDIS_ESTIMATOR_FLUX)
+ return RES_BAD_ARG;
+ SETUP_MC(flux, &estimator->fluxes[FLUX_IMPOSED]);
+ return RES_OK;
+}
+
+res_T
sdis_estimator_get_total_flux
(const struct sdis_estimator* estimator, struct sdis_mc* flux)
{
diff --git a/src/sdis_estimator_c.h b/src/sdis_estimator_c.h
@@ -31,6 +31,7 @@ enum sdis_estimator_type;
enum flux_name {
FLUX_CONVECTIVE,
FLUX_RADIATIVE,
+ FLUX_IMPOSED,
FLUX_TOTAL,
FLUX_NAMES_COUNT__
};
diff --git a/src/sdis_solve_boundary_Xd.h b/src/sdis_solve_boundary_Xd.h
@@ -446,11 +446,13 @@ XD(solve_boundary_flux)
struct accum* acc_fl = NULL; /* Per thread flux accumulator */
struct accum* acc_fc = NULL; /* Per thread convective flux accumulator */
struct accum* acc_fr = NULL; /* Per thread radiative flux accumulator */
+ struct accum* acc_fi = NULL; /* Per thread imposed flux accumulator */
size_t nrealisations = 0;
int64_t irealisation;
size_t i;
size_t view_nprims;
int progress = 0;
+ int msg1 = 0, msg2 = 0;
ATOMIC nsolved_realisations = 0;
ATOMIC res = RES_OK;
@@ -552,6 +554,7 @@ XD(solve_boundary_flux)
ALLOC_ACCUMS(acc_fc);
ALLOC_ACCUMS(acc_fl);
ALLOC_ACCUMS(acc_fr);
+ ALLOC_ACCUMS(acc_fi);
#undef ALLOC_ACCUMS
/* Create the estimator */
@@ -560,7 +563,7 @@ XD(solve_boundary_flux)
nrealisations = args->nrealisations;
omp_set_num_threads((int)scn->dev->nthreads);
- #pragma omp parallel for schedule(static)
+ #pragma omp parallel for schedule(static) shared(msg1, msg2)
for(irealisation = 0; irealisation < (int64_t)nrealisations; ++irealisation) {
struct time t0, t1;
const int ithread = omp_get_thread_num();
@@ -570,6 +573,7 @@ XD(solve_boundary_flux)
struct accum* acc_flux = &acc_fl[ithread];
struct accum* acc_fcon = &acc_fc[ithread];
struct accum* acc_frad = &acc_fr[ithread];
+ struct accum* acc_fimp = &acc_fi[ithread];
struct sXd(primitive) prim;
struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
const struct sdis_interface* interf;
@@ -578,7 +582,7 @@ XD(solve_boundary_flux)
double T_brf[3] = { 0, 0, 0 };
const double Tref = args->reference_temperature;
const double Tarad = args->ambient_radiative_temperature;
- double epsilon, hc, hr;
+ double epsilon, hc, hr, imposed_flux, imposed_temp;
size_t iprim;
double uv[DIM - 1];
float st[DIM - 1];
@@ -626,6 +630,17 @@ XD(solve_boundary_flux)
|| ( !(fmd->type == SDIS_FLUID && bmd->type == SDIS_SOLID)
&& !(fmd->type == SDIS_SOLID && bmd->type == SDIS_FLUID)))
{
+ #pragma omp critical
+ {
+ if(msg1 == 0) {
+ msg1 = 1,
+ log_err(scn->dev,
+ "%s: Attempt to compute a flux at a %s-%s interface.\n",
+ FUNC_NAME,
+ (fmd->type == SDIS_FLUID ? "fluid" : "solid"),
+ (bmd->type == SDIS_FLUID ? "fluid" : "solid"));
+ }
+ }
ATOMIC_SET(&res, RES_BAD_ARG);
continue;
}
@@ -640,14 +655,29 @@ XD(solve_boundary_flux)
/* Fetch interface parameters */
epsilon = interface_side_get_emissivity(interf, &frag);
hc = interface_get_convection_coef(interf, &frag);
-
hr = 4.0 * BOLTZMANN_CONSTANT * Tref * Tref * Tref * epsilon;
+ frag.side = solid_side;
+ imposed_flux = interface_side_get_flux(interf, &frag);
+ imposed_temp = interface_side_get_temperature(interf, &frag);
+ if(imposed_temp >= 0) {
+ /* Flux computation on T boundaries is not supported yet */
+ #pragma omp critical
+ {
+ if(msg2 == 0) {
+ msg2 = 1,
+ log_err(scn->dev,
+ "%s: Attempt to compute a flux at a Dirichlet boundary.\n",
+ FUNC_NAME);
+ }
+ }
+ ATOMIC_SET(&res, RES_BAD_ARG);
+ continue;
+ }
/* Fluid, Radiative and Solid temperatures */
flux_mask = 0;
if(hr > 0) flux_mask |= FLUX_FLAG_RADIATIVE;
if(hc > 0) flux_mask |= FLUX_FLAG_CONVECTIVE;
-
res_simul = XD(boundary_flux_realisation)(scn, rng, iprim, uv, time,
solid_side, args->fp_to_meter, Tarad, Tref, flux_mask, T_brf);
@@ -664,7 +694,8 @@ XD(solve_boundary_flux)
const double Tfluid = T_brf[2];
const double w_conv = hc * (Tboundary - Tfluid);
const double w_rad = hr * (Tboundary - Tradiative);
- const double w_total = w_conv + w_rad;
+ const double w_imp = (imposed_flux != SDIS_FLUX_NONE) ? imposed_flux : 0;
+ const double w_total = w_conv + w_rad + w_imp;
/* Temperature */
acc_temp->sum += Tboundary;
acc_temp->sum2 += Tboundary*Tboundary;
@@ -685,6 +716,10 @@ XD(solve_boundary_flux)
acc_frad->sum += w_rad;
acc_frad->sum2 += w_rad*w_rad;
++acc_frad->count;
+ /* Imposed flux */
+ acc_fimp->sum += w_imp;
+ acc_fimp->sum2 += w_imp*w_imp;
+ ++acc_fimp->count;
}
/* Update progress */
@@ -707,10 +742,12 @@ XD(solve_boundary_flux)
sum_accums(acc_fc, scn->dev->nthreads, &acc_fc[0]);
sum_accums(acc_fr, scn->dev->nthreads, &acc_fr[0]);
sum_accums(acc_fl, scn->dev->nthreads, &acc_fl[0]);
+ sum_accums(acc_fi, scn->dev->nthreads, &acc_fi[0]);
ASSERT(acc_tp[0].count == acc_fl[0].count);
ASSERT(acc_tp[0].count == acc_ti[0].count);
ASSERT(acc_tp[0].count == acc_fr[0].count);
ASSERT(acc_tp[0].count == acc_fc[0].count);
+ ASSERT(acc_tp[0].count == acc_fi[0].count);
/* Setup the estimated values */
estimator_setup_realisations_count(estimator, nrealisations, acc_tp[0].count);
@@ -718,6 +755,7 @@ XD(solve_boundary_flux)
estimator_setup_realisation_time(estimator, acc_ti[0].sum, acc_ti[0].sum2);
estimator_setup_flux(estimator, FLUX_CONVECTIVE, acc_fc[0].sum, acc_fc[0].sum2);
estimator_setup_flux(estimator, FLUX_RADIATIVE, acc_fr[0].sum, acc_fr[0].sum2);
+ estimator_setup_flux(estimator, FLUX_IMPOSED, acc_fi[0].sum, acc_fi[0].sum2);
estimator_setup_flux(estimator, FLUX_TOTAL, acc_fl[0].sum, acc_fl[0].sum2);
res = estimator_save_rng_state(estimator, rng_proxy);
@@ -733,6 +771,7 @@ exit:
if(acc_fc) MEM_RM(scn->dev->allocator, acc_fc);
if(acc_fr) MEM_RM(scn->dev->allocator, acc_fr);
if(acc_fl) MEM_RM(scn->dev->allocator, acc_fl);
+ if(acc_fi) MEM_RM(scn->dev->allocator, acc_fi);
if(scene) SXD(scene_ref_put(scene));
if(shape) SXD(shape_ref_put(shape));
if(view) SXD(scene_view_ref_put(view));
diff --git a/src/sdis_solve_probe_boundary_Xd.h b/src/sdis_solve_probe_boundary_Xd.h
@@ -346,10 +346,12 @@ XD(solve_probe_boundary_flux)
struct accum* acc_fl = NULL; /* Per thread flux accumulator */
struct accum* acc_fc = NULL; /* Per thread convective flux accumulator */
struct accum* acc_fr = NULL; /* Per thread radiative flux accumulator */
+ struct accum* acc_fi = NULL; /* Per thread imposed flux accumulator */
size_t nrealisations = 0;
int64_t irealisation = 0;
size_t i;
int progress = 0;
+ int msg1 = 0;
ATOMIC nsolved_realisations = 0;
ATOMIC res = RES_OK;
@@ -413,6 +415,11 @@ XD(solve_probe_boundary_flux)
if(!fmd || !bmd
|| ( !(fmd->type == SDIS_FLUID && bmd->type == SDIS_SOLID)
&& !(fmd->type == SDIS_SOLID && bmd->type == SDIS_FLUID))) {
+ log_err(scn->dev,
+ "%s: Attempt to compute a flux at a %s-%s interface.\n",
+ FUNC_NAME,
+ (fmd->type == SDIS_FLUID ? "fluid" : "solid"),
+ (bmd->type == SDIS_FLUID ? "fluid" : "solid"));
res = RES_BAD_ARG;
goto error;
}
@@ -449,6 +456,7 @@ XD(solve_probe_boundary_flux)
ALLOC_ACCUMS(acc_fc);
ALLOC_ACCUMS(acc_fl);
ALLOC_ACCUMS(acc_fr);
+ ALLOC_ACCUMS(acc_fi);
#undef ALLOC_ACCUMS
/* Prebuild the interface fragment */
@@ -463,7 +471,7 @@ XD(solve_probe_boundary_flux)
/* Here we go! Launch the Monte Carlo estimation */
nrealisations = args->nrealisations;
omp_set_num_threads((int)scn->dev->nthreads);
- #pragma omp parallel for schedule(static)
+ #pragma omp parallel for schedule(static) shared(msg1)
for(irealisation = 0; irealisation < (int64_t)nrealisations; ++irealisation) {
struct time t0, t1;
const int ithread = omp_get_thread_num();
@@ -473,7 +481,8 @@ XD(solve_probe_boundary_flux)
struct accum* acc_flux = &acc_fl[ithread];
struct accum* acc_fcon = &acc_fc[ithread];
struct accum* acc_frad = &acc_fr[ithread];
- double time, epsilon, hc, hr;
+ struct accum* acc_fimp = &acc_fi[ithread];
+ double time, epsilon, hc, hr, imposed_flux, imposed_temp;
int flux_mask = 0;
double T_brf[3] = { 0, 0, 0 };
const double Tref = args->reference_temperature;
@@ -494,13 +503,30 @@ XD(solve_probe_boundary_flux)
epsilon = interface_side_get_emissivity(interf, &frag);
hc = interface_get_convection_coef(interf, &frag);
hr = 4.0 * BOLTZMANN_CONSTANT * Tref * Tref * Tref * epsilon;
+ frag.side = solid_side;
+ imposed_flux = interface_side_get_flux(interf, &frag);
+ imposed_temp = interface_side_get_temperature(interf, &frag);
+ if(imposed_temp >= 0) {
+ /* Flux computation on T boundaries is not supported yet */
+ #pragma omp critical
+ {
+ if(msg1 == 0) {
+ msg1 = 1,
+ log_err(scn->dev,
+ "%s: Attempt to compute a flux at a Dirichlet boundary.\n",
+ FUNC_NAME);
+ }
+ }
+ ATOMIC_SET(&res, RES_BAD_ARG);
+ continue;
+ }
/* Fluid, Radiative and Solid temperatures */
flux_mask = 0;
if(hr > 0) flux_mask |= FLUX_FLAG_RADIATIVE;
if(hc > 0) flux_mask |= FLUX_FLAG_CONVECTIVE;
- res_simul = XD(boundary_flux_realisation)(scn, rng, args->iprim, args->uv, time,
- solid_side, args->fp_to_meter, Tarad, Tref, flux_mask, T_brf);
+ res_simul = XD(boundary_flux_realisation)(scn, rng, args->iprim, args->uv,
+ time, solid_side, args->fp_to_meter, Tarad, Tref, flux_mask, T_brf);
/* Stop time registration */
time_sub(&t0, time_current(&t1), &t0);
@@ -515,7 +541,8 @@ XD(solve_probe_boundary_flux)
const double Tfluid = T_brf[2];
const double w_conv = hc * (Tboundary - Tfluid);
const double w_rad = hr * (Tboundary - Tradiative);
- const double w_total = w_conv + w_rad;
+ const double w_imp = (imposed_flux != SDIS_FLUX_NONE) ? imposed_flux : 0;
+ const double w_total = w_conv + w_rad + w_imp;
/* Temperature */
acc_temp->sum += Tboundary;
acc_temp->sum2 += Tboundary*Tboundary;
@@ -536,6 +563,10 @@ XD(solve_probe_boundary_flux)
acc_frad->sum += w_rad;
acc_frad->sum2 += w_rad*w_rad;
++acc_frad->count;
+ /* Imposed flux */
+ acc_fimp->sum += w_imp;
+ acc_fimp->sum2 += w_imp*w_imp;
+ ++acc_fimp->count;
}
/* Update progress */
@@ -558,10 +589,12 @@ XD(solve_probe_boundary_flux)
sum_accums(acc_fc, scn->dev->nthreads, &acc_fc[0]);
sum_accums(acc_fr, scn->dev->nthreads, &acc_fr[0]);
sum_accums(acc_fl, scn->dev->nthreads, &acc_fl[0]);
+ sum_accums(acc_fi, scn->dev->nthreads, &acc_fi[0]);
ASSERT(acc_tp[0].count == acc_fl[0].count);
ASSERT(acc_tp[0].count == acc_ti[0].count);
ASSERT(acc_tp[0].count == acc_fr[0].count);
ASSERT(acc_tp[0].count == acc_fc[0].count);
+ ASSERT(acc_tp[0].count == acc_fi[0].count);
/* Setup the estimated values */
estimator_setup_realisations_count(estimator, nrealisations, acc_tp[0].count);
@@ -569,6 +602,7 @@ XD(solve_probe_boundary_flux)
estimator_setup_realisation_time(estimator, acc_ti[0].sum, acc_ti[0].sum2);
estimator_setup_flux(estimator, FLUX_CONVECTIVE, acc_fc[0].sum, acc_fc[0].sum2);
estimator_setup_flux(estimator, FLUX_RADIATIVE, acc_fr[0].sum, acc_fr[0].sum2);
+ estimator_setup_flux(estimator, FLUX_IMPOSED, acc_fi[0].sum, acc_fi[0].sum2);
estimator_setup_flux(estimator, FLUX_TOTAL, acc_fl[0].sum, acc_fl[0].sum2);
res = estimator_save_rng_state(estimator, rng_proxy);
@@ -584,6 +618,7 @@ exit:
if(acc_fc) MEM_RM(scn->dev->allocator, acc_fc);
if(acc_fr) MEM_RM(scn->dev->allocator, acc_fr);
if(acc_fl) MEM_RM(scn->dev->allocator, acc_fl);
+ if(acc_fi) MEM_RM(scn->dev->allocator, acc_fi);
if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
if(out_estimator) *out_estimator = estimator;
return (res_T)res;
