commit f794470ea50e64fb53a1312e966af01447d7b219
parent 755b4e09a5a60dd70468c2672c5657a8548916dc
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 10 Jun 2025 17:18:27 +0200
Merge branch 'feature_volrad_sw_decomposition' into develop
Diffstat:
4 files changed, 167 insertions(+), 72 deletions(-)
diff --git a/doc/htrdr-planets.1.in b/doc/htrdr-planets.1.in
@@ -20,7 +20,7 @@
.\"
.\" You should have received a copy of the GNU General Public License
.\" along with this program. If not, see <http://www.gnu.org/licenses/>.
-.Dd May 9, 2025
+.Dd May 23, 2025
.Dt HTRDR-PLANETS 1
.Os
.\""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
@@ -575,21 +575,26 @@ The other values are set to 0.
.Ss Volumic Radiative Budget
For volumic radiative budget
.Pq option Fl r
-the output is a list of 7 ASCII values per line, with as many lines as
+the output is a list of 15 ASCII values per line, with as many lines as
there are tetrahedra in the volume mesh as an argument to the
.Fl r
option.
The lines follow the order of the input meshes.
.Pp
-For each line, the first 2 values correspond to the expected value of
-the volumic radiative budget of the tetrahedron (in W/m^3) and its
-standard error.
-.Pp
-The following 3 values are the sum of the Monte Carlo weights, the sum
-of the Monte Carlo weights squared, and the total number of Monte Carlo
-weights used to calculate the preceding quantities.
-Their purpose is to help calculate the expected value and standard error
-of the volumic radiative budget for a set of tetrahedra.
+The total radiative volumic budget is decomposed into its direct and
+diffuse components.
+For each component
+.Pq total, direct and diffuse parts ,
+the following information is recorded: the average
+.Bq W/m^3 ,
+the associated standard deviation
+.Bq W/m^3 ,
+the sum of Monte Carlo weights and the sum of squared weights.
+The purpose of these last two values is to help calculate the expected value
+and the standard deviation of the volumic radiative budget for a set of
+tetrahedra.
+.Pp
+After these 12 values, the total number of realisations is recorded.
.Pp
Finally, the last two values are the estimate and associated standard
error of the calculation time per radiative path.
diff --git a/src/planets/htrdr_planets_c.h b/src/planets/htrdr_planets_c.h
@@ -68,14 +68,36 @@ struct planets_pixel_image {
HTRDR_ACCUM_NULL__ \
}
+enum planets_volrad_weight_type {
+ PLANETS_VOLRAD_TOTAL, /* 0 */
+ PLANETS_VOLRAD_DIRECT, /* 1 */
+ PLANETS_VOLRAD_DIFFUSE, /* 2 */
+ PLANETS_VOLRAD_WEIGHTS_COUNT /* 3 */
+};
+
struct planets_voxel_radiative_budget {
- struct htrdr_accum volrad_budget; /* W/m^3 */
+ struct htrdr_accum volrad_budget[PLANETS_VOLRAD_WEIGHTS_COUNT]; /* W/m^3 */
struct htrdr_accum time; /* In us */
};
-#define PLANETS_VOXEL_RADIATIVE_BUDGET { \
- HTRDR_ACCUM_NULL__, \
+#define PLANETS_VOXEL_RADIATIVE_BUDGET_NULL__ { \
+ { \
+ HTRDR_ACCUM_NULL__, \
+ HTRDR_ACCUM_NULL__, \
+ HTRDR_ACCUM_NULL__ \
+ }, \
HTRDR_ACCUM_NULL__ \
}
+static const struct planets_voxel_radiative_budget
+PLANETS_VOXEL_RADIATIVE_BUDGET_NULL = PLANETS_VOXEL_RADIATIVE_BUDGET_NULL__;
+
+enum planets_radiance_cpnt_flag {
+ PLANETS_RADIANCE_CPNT_DIRECT = BIT(0),
+ PLANETS_RADIANCE_CPNT_DIFFUSE = BIT(1),
+ PLANETS_RADIANCE_CPNT_NONE = 0,
+ PLANETS_RADIANCE_CPNT_ALL =
+ PLANETS_RADIANCE_CPNT_DIRECT
+ | PLANETS_RADIANCE_CPNT_DIFFUSE
+};
struct planets_compute_radiance_args {
struct ssp_rng* rng;
@@ -87,8 +109,11 @@ struct planets_compute_radiance_args {
double wlen; /* In nm */
size_t iband; /* Spectral band index */
size_t iquad; /* Quadrature point */
+
+ int component; /* Combination of planets_radiance_cpnt_flag */
};
-#define PLANETS_COMPUTE_RADIANCE_ARGS_NULL__ {NULL, 0, {0,0,0}, {0,0,0}, 0, 0, 0}
+#define PLANETS_COMPUTE_RADIANCE_ARGS_NULL__ \
+ {NULL, 0, {0,0,0}, {0,0,0}, 0, 0, 0, PLANETS_RADIANCE_CPNT_ALL}
static const struct planets_compute_radiance_args
PLANETS_COMPUTE_RADIANCE_ARGS_NULL = PLANETS_COMPUTE_RADIANCE_ARGS_NULL__;
diff --git a/src/planets/htrdr_planets_compute_radiance.c b/src/planets/htrdr_planets_compute_radiance.c
@@ -604,16 +604,34 @@ planets_compute_radiance
double L = 0; /* Radiance in W/m²/sr/m */
size_t nsc = 0; /* Number of surface or volume scatterings (for debug) */
int longwave = 0;
+ int shortwave = 0;
+ int direct = 0;
+ int diffuse = 0;
ASSERT(cmd && check_planets_compute_radiance_args(cmd, args) == RES_OK);
d3_set(pos, args->path_org);
d3_set(dir, args->path_dir);
+
longwave = cmd->spectral_domain.type == HTRDR_SPECTRAL_LW;
+ shortwave = !longwave;
+
+ /* In shortwave define which components are enabled */
+ if(shortwave) {
+ direct = (args->component & PLANETS_RADIANCE_CPNT_DIRECT) != 0;
+ diffuse = (args->component & PLANETS_RADIANCE_CPNT_DIFFUSE) != 0;
+ }
- if(!longwave && htrdr_planets_source_is_targeted(cmd->source, pos, dir)) {
+ /* Handle direct shortwave contribution */
+ if(shortwave
+ && direct
+ && htrdr_planets_source_is_targeted(cmd->source, pos, dir)) {
L = direct_contribution(cmd, args, pos, dir, NULL); /* In W/m²/sr/m */
}
+ /* Nothing left to do: if only the diffuse component is required
+ * in the SW, the diffuse component should not be computed */
+ if(shortwave && !diffuse) return L;
+
for(;;) {
double ev_pos[3];
double sc_dir[3];
diff --git a/src/planets/htrdr_planets_solve_volrad_budget.c b/src/planets/htrdr_planets_solve_volrad_budget.c
@@ -183,11 +183,12 @@ get_ka
return ka;
}
-static double
+static void
realisation
(struct htrdr_planets* cmd,
const struct htrdr_solve_item_args* args,
- const struct smsh_desc* volrad_mesh_desc)
+ const struct smsh_desc* volrad_mesh_desc,
+ double weights[PLANETS_VOLRAD_WEIGHTS_COUNT]) /* [W/m^3] */
{
struct planets_compute_radiance_args rad_args =
PLANETS_COMPUTE_RADIANCE_ARGS_NULL;
@@ -200,41 +201,76 @@ realisation
size_t iquad = 0; /* Quadrature point */
/* Spatial & angular integration */
+ double dir_src[3] = {0,0,0}; /* Direction toward the source */
double dir[3] = {0,0,0};
double pos[3] = {0,0,0};
+ double dir_src_pdf = 0;
+ double dir_pdf = 0;
double S = 0; /* Source [W/m^2/sr/m] */
- double L = 0; /* Radiance [W/m^2/sr/m] */
+ double L_direct = 0; /* Direct radiance [W/m^2/sr/m] */
+ double L_diffuse = 0; /* Diffuse radiance [W/m^2/sr/m] */
double ka = 0; /* Absorption coefficient */
- /* Monte Carlo weight */
- double weight = 0; /* weight [W/m^3] */
-
/* Preconditions */
ASSERT(cmd && args && volrad_mesh_desc);
ASSERT(cmd->output_type == HTRDR_PLANETS_ARGS_OUTPUT_VOLUMIC_RADIATIVE_BUDGET);
+ /* Initialise the weights */
+ memset(weights, 0, sizeof(double)*PLANETS_VOLRAD_WEIGHTS_COUNT);
+
spectral_sampling(cmd, args, &wlen, &wlen_pdf_nm, &iband, &iquad);
position_sampling(args, volrad_mesh_desc, pos);
- ssp_ran_sphere_uniform(args->rng, dir, NULL/*pdf*/);
+ ssp_ran_sphere_uniform(args->rng, dir, &dir_pdf);
+
+ S = get_source(cmd, pos, wlen); /* [W/m^2/sr/m] */
+
+ ka = get_ka(cmd, pos, iband, iquad);
+ wlen_pdf_m = wlen_pdf_nm * 1.e9; /* Transform pdf from nm^-1 to m^-1 */
/* Compute the radiance in W/m^2/sr/m */
d3_set(rad_args.path_org, pos);
- d3_set(rad_args.path_dir, dir);
rad_args.rng = args->rng;
rad_args.ithread = args->ithread;
rad_args.wlen = wlen; /* [nm] */
rad_args.iband = iband;
rad_args.iquad = iquad;
- L = planets_compute_radiance(cmd, &rad_args); /* [W/m^2/sr/m] */
- S = get_source(cmd, pos, wlen); /* [W/m^2/sr/m] */
+ if(cmd->spectral_domain.type == HTRDR_SPECTRAL_LW) {
+ /* In the longwave (radiation due to the medium), simply sample a radiative
+ * path for the sampled direction and position: the radiance is considered
+ * as purely diffuse. */
+ d3_set(rad_args.path_dir, dir);
+ L_diffuse = planets_compute_radiance(cmd, &rad_args); /* [W/m^2/sr/m] */
- ka = get_ka(cmd, pos, iband, iquad);
- wlen_pdf_m = wlen_pdf_nm * 1.e9; /* Transform pdf from nm^-1 to m^-1 */
+ /* Calculate the weights [W/m^3] */
+ weights[PLANETS_VOLRAD_DIRECT] = 0.0;
+ weights[PLANETS_VOLRAD_DIFFUSE] = ka * (L_diffuse - S) / (wlen_pdf_m * dir_pdf);
+
+ } else {
+ /* In the so-called shortwave region (actually, the radiation due the
+ * external source) is decomposed in its direct and diffuse components */
+
+ dir_src_pdf = htrdr_planets_source_sample_direction
+ (cmd->source, args->rng, pos, dir_src);
- weight = ka * (L - S) * 4*PI / wlen_pdf_m; /* [W/m^3] */
- return weight;
+ d3_set(rad_args.path_dir, dir_src);
+ rad_args.component = PLANETS_RADIANCE_CPNT_DIRECT;
+ L_direct = planets_compute_radiance(cmd, &rad_args); /* [W/m^2/sr/m] */
+
+ d3_set(rad_args.path_dir, dir);
+ rad_args.component = PLANETS_RADIANCE_CPNT_DIFFUSE;
+ L_diffuse = planets_compute_radiance(cmd, &rad_args); /* [W/m^2/sr/m] */
+
+ /* Calculate the weights [W/m^3] */
+ weights[PLANETS_VOLRAD_DIRECT] = ka * (L_direct - S) / (wlen_pdf_m * dir_src_pdf);
+ weights[PLANETS_VOLRAD_DIFFUSE] = ka * (L_diffuse - S) / (wlen_pdf_m * dir_pdf);
+ }
+
+ /* Calculate the weights [W/m^3] */
+ weights[PLANETS_VOLRAD_TOTAL] =
+ weights[PLANETS_VOLRAD_DIRECT]
+ + weights[PLANETS_VOLRAD_DIFFUSE];
}
static void
@@ -246,10 +282,6 @@ solve_volumic_radiative_budget
/* Volumic mesh on which volumic radiative budget is estimated */
struct smsh_desc volrad_mesh_desc = SMSH_DESC_NULL;
- /* Monte Carlo accumulators */
- struct htrdr_accum budget = HTRDR_ACCUM_NULL;
- struct htrdr_accum time = HTRDR_ACCUM_NULL;
-
/* Miscellaneous */
struct htrdr_planets* cmd = NULL;
struct planets_voxel_radiative_budget* voxel = data;
@@ -265,91 +297,106 @@ solve_volumic_radiative_budget
SMSH(get_desc(cmd->volrad_mesh, &volrad_mesh_desc));
+ /* Initialse voxel accumulators to 0 */
+ *voxel = PLANETS_VOXEL_RADIATIVE_BUDGET_NULL;
+
FOR_EACH(i, 0, args->nrealisations) {
/* Time recording */
struct time t0, t1;
/* Monte Carlo weights */
- double w = 0; /* [W/m^3] */
+ double w[PLANETS_VOLRAD_WEIGHTS_COUNT] = {0}; /* [W/m^3] */
double usec = 0; /* [us] */
+ int iweight = 0;
+
/* Start of realisation time recording */
time_current(&t0);
/* Run the realisation */
- w = realisation(cmd, args, &volrad_mesh_desc);
+ realisation(cmd, args, &volrad_mesh_desc, w);
/* Stop time recording */
time_sub(&t0, time_current(&t1), &t0);
usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
- /* Update the volumic radiance budget accumulator */
- budget.sum_weights += w;
- budget.sum_weights_sqr += w*w;
- budget.nweights += 1;
+ FOR_EACH(iweight, 0, PLANETS_VOLRAD_WEIGHTS_COUNT){
+ /* Update the volumic radiance budget accumulator */
+ voxel->volrad_budget[iweight].sum_weights += w[iweight];
+ voxel->volrad_budget[iweight].sum_weights_sqr += w[iweight]*w[iweight];
+ voxel->volrad_budget[iweight].nweights += 1;
+ }
/* Update the per realisation time accumulator */
- time.sum_weights += usec;
- time.sum_weights_sqr += usec*usec;
- time.nweights += 1;
+ voxel->time.sum_weights += usec;
+ voxel->time.sum_weights_sqr += usec*usec;
+ voxel->time.nweights += 1;
}
-
- /* Write voxel data */
- voxel->volrad_budget = budget;
- voxel->time = time;
}
static res_T
write_buffer
(struct htrdr_planets* cmd,
struct htrdr_buffer* buf,
- struct htrdr_accum* budget, /* W/m^3 */
- struct htrdr_accum* time, /* us */
+ struct htrdr_accum* out_budget, /* W/m^3 */
+ struct htrdr_accum* out_time, /* us */
FILE* stream)
{
struct htrdr_buffer_layout layout = HTRDR_BUFFER_LAYOUT_NULL;
size_t x = 0;
/* Preconditions */
- ASSERT(cmd && buf && budget && time && stream);
+ ASSERT(cmd && buf && out_budget && out_time && stream);
htrdr_buffer_get_layout(buf, &layout);
ASSERT(layout.height == 1);
(void)cmd;
/* Reset global accumulators */
- *budget = HTRDR_ACCUM_NULL;
- *time = HTRDR_ACCUM_NULL;
+ *out_budget = HTRDR_ACCUM_NULL;
+ *out_time = HTRDR_ACCUM_NULL;
FOR_EACH(x, 0, layout.width) {
struct planets_voxel_radiative_budget* voxel = htrdr_buffer_at(buf, x, 0);
- struct htrdr_estimate estim_budget = HTRDR_ESTIMATE_NULL;
struct htrdr_estimate estim_time = HTRDR_ESTIMATE_NULL;
-
- htrdr_accum_get_estimation(&voxel->volrad_budget, &estim_budget);
- htrdr_accum_get_estimation(&voxel->time, &estim_time);
-
- /* Write the estimate of the volumic radiative budget */
- fprintf(stream, "%g %g ", estim_budget.E, estim_budget.SE);
-
- /* Write the accumulator of the volumic radiative budget */
- fprintf(stream, "%g %g %lu ",
- voxel->volrad_budget.sum_weights,
- voxel->volrad_budget.sum_weights_sqr,
- (unsigned long)voxel->volrad_budget.nweights);
+ struct htrdr_accum* budget = NULL;
+ int iestim = 0;
+
+ budget = voxel->volrad_budget;
+ FOR_EACH(iestim, 0, PLANETS_VOLRAD_WEIGHTS_COUNT) {
+ struct htrdr_estimate estim_budget = HTRDR_ESTIMATE_NULL;
+
+ /* Write the estimate of the volumic radiative budget */
+ htrdr_accum_get_estimation(&budget[iestim], &estim_budget);
+ fprintf(stream, "%g %g ", estim_budget.E, estim_budget.SE);
+
+ /* Write the accumulator of the volumic radiative budget */
+ fprintf(stream, "%g %g ",
+ budget[iestim].sum_weights,
+ budget[iestim].sum_weights_sqr);
+ }
+
+ /* Write the number of realisations.
+ * It must be the same for all components */
+ ASSERT(budget[PLANETS_VOLRAD_TOTAL].nweights
+ == budget[PLANETS_VOLRAD_DIRECT].nweights);
+ ASSERT(budget[PLANETS_VOLRAD_TOTAL].nweights
+ == budget[PLANETS_VOLRAD_DIFFUSE].nweights);
+ fprintf(stream, "%lu ", (unsigned long)budget[PLANETS_VOLRAD_TOTAL].nweights);
/* Write the estimate of the per realisation time */
+ htrdr_accum_get_estimation(&voxel->time, &estim_time);
fprintf(stream, "%g %g\n", estim_time.E, estim_time.SE);
- /* Update the overall volumic radiative budget accumulator */
- budget->sum_weights += voxel->volrad_budget.sum_weights;
- budget->sum_weights_sqr += voxel->volrad_budget.sum_weights_sqr;
- budget->nweights += voxel->volrad_budget.nweights;
+ /* Update the overall (total) volumic radiative budget accumulator */
+ out_budget->sum_weights += budget[PLANETS_VOLRAD_TOTAL].sum_weights;
+ out_budget->sum_weights_sqr += budget[PLANETS_VOLRAD_TOTAL].sum_weights_sqr;
+ out_budget->nweights += budget[PLANETS_VOLRAD_TOTAL].nweights;
/* Update the overall per realisation time accumulator */
- time->sum_weights += voxel->time.sum_weights;
- time->sum_weights_sqr += voxel->time.sum_weights_sqr;
- time->nweights += voxel->time.nweights;
+ out_time->sum_weights += voxel->time.sum_weights;
+ out_time->sum_weights_sqr += voxel->time.sum_weights_sqr;
+ out_time->nweights += voxel->time.nweights;
}
return RES_OK;
}
