stardis-solver

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

test_sdis_solve_medium.c (19082B)


      1 /* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com)
      2  *
      3  * This program is free software: you can redistribute it and/or modify
      4  * it under the terms of the GNU General Public License as published by
      5  * the Free Software Foundation, either version 3 of the License, or
      6  * (at your option) any later version.
      7  *
      8  * This program is distributed in the hope that it will be useful,
      9  * but WITHOUT ANY WARRANTY; without even the implied warranty of
     10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
     11  * GNU General Public License for more details.
     12  *
     13  * You should have received a copy of the GNU General Public License
     14  * along with this program. If not, see <http://www.gnu.org/licenses/>. */
     15 
     16 #include "sdis.h"
     17 #include "test_sdis_utils.h"
     18 
     19 #include <rsys/math.h>
     20 #include <rsys/stretchy_array.h>
     21 #include <star/ssp.h>
     22 #include <star/s3dut.h>
     23 
     24 #include <string.h>
     25 
     26 #define Tf0 300.0
     27 #define Tf1 330.0
     28 #define N 1000ul /* #realisations */
     29 #define Np 10000ul /* #realisations precise */
     30 
     31 /*
     32  * The scene is composed of 2 super shapes whose temperature is unknown. The
     33  * first super shape is surrounded by a fluid whose temperature is Tf0 while
     34  * the second one is in fluid whose temperature is Tf1. The temperatures of the
     35  * super shape 0 and 1 are thus uniform and equal to Tf0 and Tf1, respectively.
     36  *
     37  * This program performs 2 tests. In the first one, the super shapes 0 and 1
     38  * have different media; the medium solver thus estimates the
     39  * temperature of one super shape. In the second test, the scene is updated to
     40  * use the same medium for the 2 super shapes. In this case, when invoked on
     41  * the right medium, the estimated temperature T is equal to :
     42  *
     43  *    T = Tf0 * V0/(V0 + V1) + Tf1 * V1/(V0 + V1)
     44  *
     45  * with V0 and V1 the volume of the super shapes 0 and 1, respectively.
     46  */
     47 
     48 /*******************************************************************************
     49  * Geometry
     50  ******************************************************************************/
     51 struct context {
     52   struct s3dut_mesh_data msh0;
     53   struct s3dut_mesh_data msh1;
     54   struct sdis_interface* interf0;
     55   struct sdis_interface* interf1;
     56 };
     57 
     58 static void
     59 get_indices(const size_t itri, size_t ids[3], void* context)
     60 {
     61   const struct context* ctx = context;
     62   /* Note that we swap the indices to ensure that the triangle normals point
     63    * inward the super shape */
     64   if(itri < ctx->msh0.nprimitives) {
     65     ids[0] = ctx->msh0.indices[itri*3+0];
     66     ids[2] = ctx->msh0.indices[itri*3+1];
     67     ids[1] = ctx->msh0.indices[itri*3+2];
     68   } else {
     69     const size_t itri2 = itri - ctx->msh0.nprimitives;
     70     ids[0] = ctx->msh1.indices[itri2*3+0] + ctx->msh0.nvertices;
     71     ids[2] = ctx->msh1.indices[itri2*3+1] + ctx->msh0.nvertices;
     72     ids[1] = ctx->msh1.indices[itri2*3+2] + ctx->msh0.nvertices;
     73   }
     74 }
     75 
     76 static void
     77 get_position(const size_t ivert, double pos[3], void* context)
     78 {
     79   const struct context* ctx = context;
     80   if(ivert < ctx->msh0.nvertices) {
     81     pos[0] = ctx->msh0.positions[ivert*3+0] - 2.0;
     82     pos[1] = ctx->msh0.positions[ivert*3+1];
     83     pos[2] = ctx->msh0.positions[ivert*3+2];
     84   } else {
     85     const size_t ivert2 = ivert - ctx->msh0.nvertices;
     86     pos[0] = ctx->msh1.positions[ivert2*3+0] + 2.0;
     87     pos[1] = ctx->msh1.positions[ivert2*3+1];
     88     pos[2] = ctx->msh1.positions[ivert2*3+2];
     89   }
     90 }
     91 
     92 static void
     93 get_interface(const size_t itri, struct sdis_interface** bound, void* context)
     94 {
     95   const struct context* ctx = context;
     96   *bound = itri < ctx->msh0.nprimitives ? ctx->interf0 : ctx->interf1;
     97 }
     98 
     99 /*******************************************************************************
    100  * Fluid medium
    101  ******************************************************************************/
    102 struct fluid {
    103   double temperature;
    104 };
    105 
    106 static double
    107 fluid_get_temperature
    108   (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
    109 {
    110   CHK(data != NULL && vtx != NULL);
    111   return ((const struct fluid*)sdis_data_cget(data))->temperature;
    112 }
    113 
    114 /*******************************************************************************
    115  * Solid medium
    116  ******************************************************************************/
    117 struct solid {
    118   double cp;
    119   double lambda;
    120   double rho;
    121   double delta;
    122   double temperature;
    123 };
    124 
    125 static double
    126 solid_get_calorific_capacity
    127   (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
    128 {
    129   CHK(data != NULL && vtx != NULL);
    130   return ((const struct solid*)sdis_data_cget(data))->cp;
    131 }
    132 
    133 static double
    134 solid_get_thermal_conductivity
    135   (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
    136 {
    137   CHK(data != NULL && vtx != NULL);
    138   return ((const struct solid*)sdis_data_cget(data))->lambda;
    139 }
    140 
    141 static double
    142 solid_get_volumic_mass
    143   (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
    144 {
    145   CHK(data != NULL && vtx != NULL);
    146   return ((const struct solid*)sdis_data_cget(data))->rho;
    147 }
    148 
    149 static double
    150 solid_get_delta
    151   (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
    152 {
    153   CHK(data != NULL && vtx != NULL);
    154   return ((const struct solid*)sdis_data_cget(data))->delta;
    155 }
    156 
    157 static double
    158 solid_get_temperature
    159   (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
    160 {
    161   CHK(data != NULL && vtx != NULL);
    162   return ((const struct solid*)sdis_data_cget(data))->temperature;
    163 }
    164 
    165 struct interf {
    166   double hc;
    167   double epsilon;
    168   double specular_fraction;
    169 };
    170 
    171 static double
    172 interface_get_convection_coef
    173   (const struct sdis_interface_fragment* frag, struct sdis_data* data)
    174 {
    175   CHK(data != NULL && frag != NULL);
    176   return ((const struct interf*)sdis_data_cget(data))->hc;
    177 }
    178 
    179 static double
    180 interface_get_emissivity
    181   (const struct sdis_interface_fragment* frag,
    182    const unsigned source_id,
    183    struct sdis_data* data)
    184 {
    185   (void)source_id;
    186   CHK(data != NULL && frag != NULL);
    187   return ((const struct interf*)sdis_data_cget(data))->epsilon;
    188 }
    189 
    190 static double
    191 interface_get_specular_fraction
    192   (const struct sdis_interface_fragment* frag,
    193    const unsigned source_id,
    194    struct sdis_data* data)
    195 {
    196   (void)source_id;
    197   CHK(data != NULL && frag != NULL);
    198   return ((const struct interf*)sdis_data_cget(data))->specular_fraction;
    199 }
    200 
    201 /*******************************************************************************
    202  * Test
    203  ******************************************************************************/
    204 int
    205 main(int argc, char** argv)
    206 {
    207   struct s3dut_super_formula f0 = S3DUT_SUPER_FORMULA_NULL;
    208   struct s3dut_super_formula f1 = S3DUT_SUPER_FORMULA_NULL;
    209   struct s3dut_mesh* msh0 = NULL;
    210   struct s3dut_mesh* msh1 = NULL;
    211   struct sdis_mc T = SDIS_MC_NULL;
    212   struct sdis_mc T2 = SDIS_MC_NULL;
    213   struct sdis_mc time = SDIS_MC_NULL;
    214   struct sdis_device* dev = NULL;
    215   struct sdis_medium* solid0 = NULL;
    216   struct sdis_medium* solid1 = NULL;
    217   struct sdis_medium* fluid0 = NULL;
    218   struct sdis_medium* fluid1 = NULL;
    219   struct sdis_interface* solid0_fluid0 = NULL;
    220   struct sdis_interface* solid0_fluid1 = NULL;
    221   struct sdis_interface* solid1_fluid1 = NULL;
    222   struct sdis_scene* scn = NULL;
    223   struct sdis_data* data = NULL;
    224   struct sdis_estimator* estimator = NULL;
    225   struct sdis_estimator* estimator2 = NULL;
    226   struct sdis_green_function* green = NULL;
    227   struct fluid* fluid_param = NULL;
    228   struct solid* solid_param = NULL;
    229   struct interf* interface_param = NULL;
    230   struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
    231   struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
    232   struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
    233   struct sdis_interface_shader interface_shader = SDIS_INTERFACE_SHADER_NULL;
    234   struct sdis_solve_medium_args solve_args = SDIS_SOLVE_MEDIUM_ARGS_DEFAULT;
    235   struct ssp_rng* rng = NULL;
    236   struct context ctx;
    237   double ref = 0;
    238   double v, v0, v1;
    239   size_t nreals;
    240   size_t nfails;
    241   size_t ntris;
    242   size_t nverts;
    243   int is_master_process;
    244   (void)argc, (void)argv;
    245 
    246   create_default_device(&argc, &argv, &is_master_process, &dev);
    247 
    248   fluid_shader.temperature = fluid_get_temperature;
    249 
    250   /* Create the fluid0 medium */
    251   OK(sdis_data_create
    252     (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data));
    253   fluid_param = sdis_data_get(data);
    254   fluid_param->temperature = Tf0;
    255   OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid0));
    256   OK(sdis_data_ref_put(data));
    257 
    258   /* Create the fluid1 medium */
    259   OK(sdis_data_create
    260     (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data));
    261   fluid_param = sdis_data_get(data);
    262   fluid_param->temperature = Tf1;
    263   OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid1));
    264   OK(sdis_data_ref_put(data));
    265 
    266   /* Setup the solid shader */
    267   solid_shader.calorific_capacity = solid_get_calorific_capacity;
    268   solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
    269   solid_shader.volumic_mass = solid_get_volumic_mass;
    270   solid_shader.delta = solid_get_delta;
    271   solid_shader.temperature = solid_get_temperature;
    272 
    273   /* Create the solid0 medium */
    274   OK(sdis_data_create
    275     (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data));
    276   solid_param = sdis_data_get(data);
    277   solid_param->cp = 1.0;
    278   solid_param->lambda = 0.1;
    279   solid_param->rho = 1.0;
    280   solid_param->delta = 1.0/20.0;
    281   solid_param->temperature = SDIS_TEMPERATURE_NONE; /* Unknown temperature */
    282   OK(sdis_solid_create(dev, &solid_shader, data, &solid0));
    283   OK(sdis_data_ref_put(data));
    284 
    285   /* Create the solid1 medium */
    286   OK(sdis_data_create
    287     (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data));
    288   solid_param = sdis_data_get(data);
    289   solid_param->cp = 1.0;
    290   solid_param->lambda = 1.0;
    291   solid_param->rho = 1.0;
    292   solid_param->delta = 1.0/20.0;
    293   solid_param->temperature = SDIS_TEMPERATURE_NONE; /* Unknown temperature */
    294   OK(sdis_solid_create(dev, &solid_shader, data, &solid1));
    295   OK(sdis_data_ref_put(data));
    296 
    297   /* Create the interfaces */
    298   OK(sdis_data_create(dev, sizeof(struct interf),
    299     ALIGNOF(struct interf), NULL, &data));
    300   interface_param = sdis_data_get(data);
    301   interface_param->hc = 0.5;
    302   interface_param->epsilon = 0;
    303   interface_param->specular_fraction = 0;
    304   interface_shader.convection_coef = interface_get_convection_coef;
    305   interface_shader.front = SDIS_INTERFACE_SIDE_SHADER_NULL;
    306   interface_shader.back.temperature = NULL;
    307   interface_shader.back.emissivity = interface_get_emissivity;
    308   interface_shader.back.specular_fraction = interface_get_specular_fraction;
    309   OK(sdis_interface_create
    310     (dev, solid0, fluid0, &interface_shader, data, &solid0_fluid0));
    311   OK(sdis_interface_create
    312     (dev, solid0, fluid1, &interface_shader, data, &solid0_fluid1));
    313   OK(sdis_interface_create
    314     (dev, solid1, fluid1, &interface_shader, data, &solid1_fluid1));
    315   OK(sdis_data_ref_put(data));
    316 
    317   /* Create the mesh0 */
    318   f0.A = 1; f0.B = 1; f0.M = 3; f0.N0 = 1; f0.N1 = 1; f0.N2 = 2;
    319   f1.A = 1; f1.B = 1; f1.M = 10; f1.N0 = 1; f1.N1 = 1; f1.N2 = 3;
    320   OK(s3dut_create_super_shape(NULL, &f0, &f1, 1, 64, 32, &msh0));
    321   OK(s3dut_mesh_get_data(msh0, &ctx.msh0));
    322 
    323   /* Create the mesh1 */
    324   f0.A = 1; f0.B = 1; f0.M = 10; f0.N0 = 1; f0.N1 = 1; f0.N2 = 5;
    325   f1.A = 1; f1.B = 1; f1.M = 1; f1.N0 = 1; f1.N1 = 1; f1.N2 = 1;
    326   OK(s3dut_create_super_shape(NULL, &f0, &f1, 1, 64, 32, &msh1));
    327   OK(s3dut_mesh_get_data(msh1, &ctx.msh1));
    328 
    329   /* Create the scene */
    330   ctx.interf0 = solid0_fluid0;
    331   ctx.interf1 = solid1_fluid1;
    332   ntris = ctx.msh0.nprimitives + ctx.msh1.nprimitives;
    333   nverts = ctx.msh0.nvertices + ctx.msh1.nvertices;
    334 #if 0
    335   {
    336     double* vertices = NULL;
    337     size_t* indices = NULL;
    338     size_t i;
    339     CHK(vertices = MEM_CALLOC(&allocator, nverts*3, sizeof(*vertices)));
    340     CHK(indices = MEM_CALLOC(&allocator, ntris*3, sizeof(*indices)));
    341     FOR_EACH(i, 0, ntris) get_indices(i, indices + i*3, &ctx);
    342     FOR_EACH(i, 0, nverts) get_position(i, vertices + i*3, &ctx);
    343     dump_mesh(stdout, vertices, nverts, indices, ntris);
    344     MEM_RM(&allocator, vertices);
    345     MEM_RM(&allocator, indices);
    346   }
    347 #endif
    348 
    349   scn_args.get_indices = get_indices;
    350   scn_args.get_interface = get_interface;
    351   scn_args.get_position = get_position;
    352   scn_args.nprimitives = ntris;
    353   scn_args.nvertices = nverts;
    354   scn_args.context = &ctx;
    355   OK(sdis_scene_create(dev, &scn_args, &scn));
    356 
    357   BA(sdis_scene_get_medium_spread(NULL, solid0, &v0));
    358   BA(sdis_scene_get_medium_spread(scn, NULL, &v0));
    359   BA(sdis_scene_get_medium_spread(scn, solid0, NULL));
    360   OK(sdis_scene_get_medium_spread(scn, solid0, &v0));
    361   CHK(v0 > 0);
    362   OK(sdis_scene_get_medium_spread(scn, solid1, &v1));
    363   CHK(v1 > 0);
    364   OK(sdis_scene_get_medium_spread(scn, fluid0, &v));
    365   CHK(v == 0);
    366   OK(sdis_scene_get_medium_spread(scn, fluid1, &v));
    367   CHK(v == 0);
    368 
    369   solve_args.nrealisations = N;
    370   solve_args.medium = solid0;
    371   solve_args.time_range[0] = INF;
    372   solve_args.time_range[1] = INF;
    373 
    374   BA(sdis_solve_medium(NULL, &solve_args, &estimator));
    375   BA(sdis_solve_medium(scn, NULL, &estimator));
    376   BA(sdis_solve_medium(scn, &solve_args, NULL));
    377   solve_args.nrealisations = 0;
    378   BA(sdis_solve_medium(scn, &solve_args, &estimator));
    379   solve_args.nrealisations = N;
    380   solve_args.medium = NULL;
    381   BA(sdis_solve_medium(scn, &solve_args, &estimator));
    382   solve_args.medium = solid0;
    383   solve_args.time_range[0] = solve_args.time_range[1] = -1;
    384   BA(sdis_solve_medium(scn, &solve_args, &estimator));
    385   solve_args.time_range[0] = 1;
    386   BA(sdis_solve_medium(scn, &solve_args, &estimator));
    387   solve_args.time_range[1] = 0;
    388   BA(sdis_solve_medium(scn, &solve_args, &estimator));
    389   solve_args.time_range[0] = solve_args.time_range[1] = INF;
    390   OK(sdis_solve_medium(scn, &solve_args, &estimator));
    391 
    392   if(!is_master_process) {
    393     CHK(estimator == NULL);
    394   } else {
    395     OK(sdis_estimator_get_realisation_count(estimator, &nreals));
    396     OK(sdis_estimator_get_failure_count(estimator, &nfails));
    397     OK(sdis_estimator_get_temperature(estimator, &T));
    398     OK(sdis_estimator_get_realisation_time(estimator, &time));
    399     printf("Shape0 temperature = "STR(Tf0)" ~ %g +/- %g\n", T.E, T.SE);
    400     printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE);
    401     printf("#failures = %lu/%lu\n\n", (unsigned long)nfails, N);
    402     CHK(eq_eps(T.E, Tf0, T.SE));
    403     CHK(nreals + nfails == N);
    404     OK(sdis_estimator_ref_put(estimator));
    405   }
    406 
    407   solve_args.medium = solid1;
    408 
    409   /* Check simulation error handling when paths are registered */
    410   solve_args.nrealisations = 10;
    411   solve_args.register_paths = SDIS_HEAT_PATH_ALL;
    412   fluid_param->temperature = SDIS_TEMPERATURE_NONE;
    413   BA(sdis_solve_medium(scn, &solve_args, &estimator));
    414   fluid_param->temperature = Tf1;
    415   OK(sdis_solve_medium(scn, &solve_args, &estimator));
    416   if(is_master_process) {
    417     OK(sdis_estimator_ref_put(estimator));
    418   }
    419   solve_args.nrealisations = N;
    420   solve_args.register_paths = SDIS_HEAT_PATH_NONE;
    421 
    422   OK(sdis_solve_medium(scn, &solve_args, &estimator));
    423   if(is_master_process) {
    424     OK(sdis_estimator_get_realisation_count(estimator, &nreals));
    425     OK(sdis_estimator_get_failure_count(estimator, &nfails));
    426     OK(sdis_estimator_get_temperature(estimator, &T));
    427     OK(sdis_estimator_get_realisation_time(estimator, &time));
    428     printf("Shape1 temperature = "STR(Tf1)" ~ %g +/- %g\n", T.E, T.SE);
    429     printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE);
    430     printf("#failures = %lu/%lu\n\n", (unsigned long)nfails, N);
    431     CHK(eq_eps(T.E, Tf1, T.SE));
    432     CHK(nreals + nfails == N);
    433     OK(sdis_estimator_ref_put(estimator));
    434   }
    435 
    436   /* Create a new scene with the same medium in the 2 super shapes */
    437   OK(sdis_scene_ref_put(scn));
    438   ctx.interf0 = solid0_fluid0;
    439   ctx.interf1 = solid0_fluid1;
    440   OK(sdis_scene_create(dev, &scn_args, &scn));
    441 
    442   OK(sdis_scene_get_medium_spread(scn, solid0, &v));
    443   CHK(eq_eps(v, v0+v1, 1.e-6));
    444 
    445   solve_args.medium = solid1;
    446   BA(sdis_solve_medium(scn, &solve_args, &estimator));
    447   solve_args.medium = solid0;
    448   solve_args.nrealisations = Np;
    449   OK(sdis_solve_medium(scn, &solve_args, &estimator));
    450   if(is_master_process) {
    451     OK(sdis_estimator_get_temperature(estimator, &T));
    452     OK(sdis_estimator_get_realisation_time(estimator, &time));
    453     OK(sdis_estimator_get_realisation_count(estimator, &nreals));
    454     OK(sdis_estimator_get_failure_count(estimator, &nfails));
    455     ref = Tf0 * v0/v + Tf1 * v1/v;
    456     printf("Shape0 + Shape1 temperature = %g ~ %g +/- %g\n", ref, T.E, T.SE);
    457     printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE);
    458     printf("#failures = %lu/%lu\n", (unsigned long)nfails, Np);
    459     CHK(eq_eps(T.E, ref, T.SE*3));
    460   }
    461 
    462   /* Check RNG type */
    463   solve_args.rng_state = NULL;
    464   solve_args.rng_type = SSP_RNG_TYPE_NULL;
    465   BA(sdis_solve_medium(scn, &solve_args, &estimator2));
    466   solve_args.rng_type =
    467     SDIS_SOLVE_MEDIUM_ARGS_DEFAULT.rng_type == SSP_RNG_THREEFRY
    468     ? SSP_RNG_MT19937_64 : SSP_RNG_THREEFRY;
    469   OK(sdis_solve_medium(scn, &solve_args, &estimator2));
    470   if(is_master_process) {
    471     OK(sdis_estimator_get_temperature(estimator2, &T2));
    472     CHK(eq_eps(T2.E, ref, 3*T2.SE));
    473     CHK(T2.E != T.E);
    474     OK(sdis_estimator_ref_put(estimator2));
    475   }
    476 
    477   /* Check RNG state */
    478   OK(ssp_rng_create(NULL, SSP_RNG_THREEFRY, &rng));
    479   OK(ssp_rng_discard(rng, 31415926535)); /* Move the RNG state  */
    480   solve_args.rng_state = rng;
    481   solve_args.rng_type = SSP_RNG_TYPE_NULL;
    482   OK(sdis_solve_medium(scn, &solve_args, &estimator2));
    483   OK(ssp_rng_ref_put(rng));
    484   if(is_master_process) {
    485     OK(sdis_estimator_get_temperature(estimator2, &T2));
    486     CHK(eq_eps(T2.E, ref, 3*T2.SE));
    487     CHK(T2.E != T.E);
    488     OK(sdis_estimator_ref_put(estimator2));
    489   }
    490 
    491   /* Restore args */
    492   solve_args.rng_state = SDIS_SOLVE_PROBE_ARGS_DEFAULT.rng_state;
    493   solve_args.rng_type = SDIS_SOLVE_PROBE_ARGS_DEFAULT.rng_type;
    494 
    495 
    496 
    497   /* Solve green */
    498   BA(sdis_solve_medium_green_function(NULL, &solve_args, &green));
    499   BA(sdis_solve_medium_green_function(scn, NULL, &green));
    500   BA(sdis_solve_medium_green_function(scn, &solve_args, NULL));
    501   solve_args.nrealisations = 0;
    502   BA(sdis_solve_medium_green_function(scn, &solve_args, &green));
    503   solve_args.nrealisations = Np;
    504   solve_args.medium = NULL;
    505   BA(sdis_solve_medium_green_function(scn, &solve_args, &green));
    506   solve_args.medium = solid1;
    507   BA(sdis_solve_medium_green_function(scn, &solve_args, &green));
    508   solve_args.medium = solid0;
    509   solve_args.picard_order = 0;
    510   BA(sdis_solve_medium_green_function(scn, &solve_args, &green));
    511   solve_args.picard_order = 1;
    512   OK(sdis_solve_medium_green_function(scn, &solve_args, &green));
    513 
    514   if(!is_master_process) {
    515     CHK(green == NULL);
    516   } else {
    517     OK(sdis_green_function_solve(green, &estimator2));
    518     check_green_function(green);
    519     check_estimator_eq(estimator, estimator2);
    520     check_green_serialization(green, scn);
    521 
    522     OK(sdis_green_function_ref_put(green));
    523 
    524     OK(sdis_estimator_ref_put(estimator));
    525     OK(sdis_estimator_ref_put(estimator2));
    526   }
    527 
    528   /* Release */
    529   OK(s3dut_mesh_ref_put(msh0));
    530   OK(s3dut_mesh_ref_put(msh1));
    531   OK(sdis_medium_ref_put(fluid0));
    532   OK(sdis_medium_ref_put(fluid1));
    533   OK(sdis_medium_ref_put(solid0));
    534   OK(sdis_medium_ref_put(solid1));
    535   OK(sdis_interface_ref_put(solid0_fluid0));
    536   OK(sdis_interface_ref_put(solid0_fluid1));
    537   OK(sdis_interface_ref_put(solid1_fluid1));
    538   OK(sdis_scene_ref_put(scn));
    539   free_default_device(dev);
    540 
    541   CHK(mem_allocated_size() == 0);
    542   return 0;
    543 }