stardis-solver

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

commit 1ab9726c0f55a151756d1b0512dce35dc1d32e24
parent 3dd4f89f76e03db1779d5460ebb41fe1c5c583de
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Fri, 12 Apr 2019 15:56:23 +0200

Improve the accuracy of the conductive random walk

Approximate the distance to the nearest boundary by tracing 4 or 6 rays
whether the scene is 2D or 3D, respectively.

Currently this new conductive random walk does not properly handle the
corrective volumetric power term.

Diffstat:
Mcmake/CMakeLists.txt | 2++
Msrc/sdis_heat_path_conductive_Xd.h | 169+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------
Msrc/sdis_realisation_Xd.h | 4++--
Msrc/sdis_scene_Xd.h | 6+++---
Asrc/test_sdis_solid_random_walk_robustness.c | 333+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Msrc/test_sdis_solve_probe3.c | 4++--
Msrc/test_sdis_volumic_power4_2d.c | 2+-
7 files changed, 478 insertions(+), 42 deletions(-)

diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt @@ -181,6 +181,7 @@ if(NOT NO_TEST) new_test(test_sdis_interface) new_test(test_sdis_medium) new_test(test_sdis_scene) + new_test(test_sdis_solid_random_walk_robustness) new_test(test_sdis_solve_camera) new_test(test_sdis_solve_probe) new_test(test_sdis_solve_probe2) @@ -206,6 +207,7 @@ if(NOT NO_TEST) add_test(test_sdis_volumic_power3_2d test_sdis_volumic_power3_2d) endif() + target_link_libraries(test_sdis_solid_random_walk_robustness Star3DUT) target_link_libraries(test_sdis_solve_medium Star3DUT) target_link_libraries(test_sdis_solve_probe3 Star3DUT) target_link_libraries(test_sdis_solve_probe3_2d ${MATH_LIB}) diff --git a/src/sdis_heat_path_conductive_Xd.h b/src/sdis_heat_path_conductive_Xd.h @@ -24,6 +24,130 @@ #include "sdis_Xd_begin.h" +/******************************************************************************* + * Helper functions + ******************************************************************************/ +/* Sample the next direction to walk toward and compute the distance to travel. + * Return the sampled direction `dir0', the distance to travel along this + * direction, the hit `hit0' along `dir0' wrt to the returned distance, the + * direction `dir1' used to adjust the displacement distance, and the hit + * `hit1' along `dir1' used to adjust the displacement distance. */ +static float +XD(sample_next_step) + (struct sdis_scene* scn, + struct ssp_rng* rng, + const float pos[DIM], + const float delta_solid, + float dir0[DIM], /* Sampled direction */ + float dir1[DIM], /* Direction used to adjust delta */ + struct sXd(hit)* hit0, /* Hit along the sampled direction */ + struct sXd(hit)* hit1) /* Hit used to adjust delta */ +{ + float dirs[2*DIM][DIM]; + float range[2]; + float delta; + struct sXd(hit) hit= SXD_HIT_NULL; + int idir; + int idir1; + ASSERT(scn && rng && pos && delta_solid>0 && dir0 && dir1 && hit0 && hit1); + + *hit0 = SXD_HIT_NULL; + *hit1 = SXD_HIT_NULL; + +#if DIM == 2 + /* Sample a main direction around 2PI */ + ssp_ran_circle_uniform_float(rng, dirs[0], NULL); + + /* Compute in dirs[2] a direction orthogonal to dirs[0] */ + dirs[2][0] = -dirs[0][1]; + dirs[2][1] = dirs[0][0]; + ASSERT(f2_is_normalized(dirs[2])); + ASSERT(eq_epsf(f2_dot(dirs[0], dirs[2]), 0, 1.e-6f)); + + /* Negate the orthornormal frame */ + f2_minus(dirs[1], dirs[0]); + f2_minus(dirs[3], dirs[2]); +#else + { + float dir_abs[DIM]; + int i, j, k; + + /* Sample a main direction around 4PI */ + ssp_ran_sphere_uniform_float(rng, dirs[0], NULL); + + /* Find the index of the maximum coordinate of the sampled direction */ + dir_abs[0] = absf(dirs[0][0]); + dir_abs[1] = absf(dirs[0][1]); + dir_abs[2] = absf(dirs[0][2]); + i = dir_abs[0] > dir_abs[1] + ? (dir_abs[0] > dir_abs[2] ? 0 : 2) + : (dir_abs[1] > dir_abs[2] ? 1 : 2); + j = (i+1) % 3; + k = (j+1) % 3; + + /* Compute a direction orthogonal to the sample dir */ + dirs[2][i] = -(dirs[0][j]*dirs[0][j] + dirs[0][k]*dirs[0][k]) / dirs[0][i]; + dirs[2][j] = dirs[0][j]; + dirs[2][k] = dirs[0][k]; + f3_normalize(dirs[2], dirs[2]); + + /* Complete the orthonormal frame */ + f3_cross(dirs[4], dirs[0], dirs[2]); + f3_normalize(dirs[4], dirs[4]); + + /* Negate the orthonormal frame */ + f3_minus(dirs[1], dirs[0]); + f3_minus(dirs[3], dirs[2]); + f3_minus(dirs[5], dirs[4]); + + ASSERT(f3_is_normalized(dirs[2])); + ASSERT(f3_is_normalized(dirs[4])); + ASSERT(eq_epsf(f3_dot(dirs[0], dirs[2]), 0, 1.e-6f)); + ASSERT(eq_epsf(f3_dot(dirs[0], dirs[4]), 0, 1.e-6f)); + ASSERT(eq_epsf(f3_dot(dirs[2], dirs[4]), 0, 1.e-6f)); + } +#endif + + /* Use the previously computed orthornormal frame to estimate the minimum + * distance from `pos' to the scene boundary */ + range[0] = 0.f; + range[1] = delta_solid*RAY_RANGE_MAX_SCALE; + delta = FLT_MAX; + idir1 = 0; + FOR_EACH(idir, 0, 2*DIM) { + SXD(scene_view_trace_ray(scn->sXd(view), pos, dirs[idir], range, NULL, &hit)); + if(idir == 0) *hit0 = hit; + if(hit.distance < delta) { + delta = hit.distance; + *hit1 = hit; + idir1 = idir; + } + } + + if(delta == FLT_MAX) { + /* Hit nothing along all tested directions. Set delta to delta_solid. */ + delta = delta_solid; + } else if + ( delta != hit0->distance + && eq_eps(hit0->distance, delta, delta_solid*(RAY_RANGE_MAX_SCALE-1))) { + /* Set delta to the main hit distance if it is roughly equal to it in order + * to avoid numerical issues on moving along the main direction. Use the + * RAY_RANGE_MAX_SCALE factor to define the `epsilon' used by this + * comparison */ + delta = hit0->distance; + *hit1 = *hit0; + idir1 = 0; + } + + fX(set)(dir0, dirs[0]); + fX(set)(dir1, dirs[idir1]); + + return delta; +} + +/******************************************************************************* + * Local function + ******************************************************************************/ res_T XD(conductive_path) (struct sdis_scene* scn, @@ -70,7 +194,6 @@ XD(conductive_path) double power_factor = 0; double power; float delta, delta_solid; /* Random walk numerical parameter */ - float range[2]; float dir0[DIM], dir1[DIM]; float org[DIM]; @@ -109,39 +232,19 @@ XD(conductive_path) goto error; } -#if DIM == 2 - /* Sample a direction around 2PI */ - ssp_ran_circle_uniform_float(rng, dir0, NULL); -#else - /* Sample a direction around 4PI */ - ssp_ran_sphere_uniform_float(rng, dir0, NULL); -#endif - - /* Trace a ray along the sampled direction and its opposite to check if a - * surface is hit in [0, delta_solid]. */ fX_set_dX(org, rwalk->vtx.P); - fX(minus)(dir1, dir0); - hit0 = hit1 = SXD_HIT_NULL; - range[0] = 0.f, range[1] = delta_solid*RAY_RANGE_MAX_SCALE; - SXD(scene_view_trace_ray(scn->sXd(view), org, dir0, range, NULL, &hit0)); - SXD(scene_view_trace_ray(scn->sXd(view), org, dir1, range, NULL, &hit1)); - - if(SXD_HIT_NONE(&hit0) && SXD_HIT_NONE(&hit1)) { - /* Hit nothing: move along dir0 of the original delta */ - delta = delta_solid; - - /* Add the volumic power density to the measured temperature */ - if(power != SDIS_VOLUMIC_POWER_NONE) { + + /* Sample the direction to walk toward and compute the distance to travel */ + delta = XD(sample_next_step) + (scn, rng, org, delta_solid, dir0, dir1, &hit0, &hit1); + + /* Add the volumic power density to the measured temperature */ + if(power != SDIS_VOLUMIC_POWER_NONE) { + if((S3D_HIT_NONE(&hit0) && S3D_HIT_NONE(&hit1))) { /* Hit nothing */ const double delta_in_meter = delta * fp_to_meter; power_factor = delta_in_meter * delta_in_meter / (2.0 * DIM * lambda); T->value += power * power_factor; - } - } else { - /* Hit something: move along dir0 of the minimum hit distance */ - delta = MMIN(hit0.distance, hit1.distance); - - /* Add the volumic power density to the measured temperature */ - if(power != SDIS_VOLUMIC_POWER_NONE) { + } else { const double delta_s_adjusted = delta_solid * RAY_RANGE_MAX_SCALE; const double delta_s_in_meter = delta_solid * fp_to_meter; double h; @@ -226,10 +329,8 @@ XD(conductive_path) } } - /* Define if the random walk hits something along dir0. Multiply delta by - * the empirical ray range scale factor to ensure that once moved, the - * random walk does not lie in the uncertainty zone near the geometry */ - if(hit0.distance > delta * RAY_RANGE_MAX_SCALE) { + /* Define if the random walk hits something along dir0 */ + if(hit0.distance > delta) { rwalk->hit = SXD_HIT_NULL; rwalk->hit_side = SDIS_SIDE_NULL__; } else { diff --git a/src/sdis_realisation_Xd.h b/src/sdis_realisation_Xd.h @@ -45,9 +45,9 @@ XD(compute_temperature) }* stack = NULL; size_t istack = 0; #endif - /* Maximum accepted #failures before stopping the realisation */ struct sdis_heat_vertex* heat_vtx = NULL; - const size_t MAX_FAILS = 10; + /* Maximum accepted #failures before stopping the realisation */ + const size_t MAX_FAILS = 1; res_T res = RES_OK; ASSERT(scn && fp_to_meter > 0 && ctx && rwalk && rng && T); diff --git a/src/sdis_scene_Xd.h b/src/sdis_scene_Xd.h @@ -850,8 +850,8 @@ XD(scene_get_medium) goto error; } } - /* Discard the hit if it is on a vertex, i.e. between 2 segments, and - * target a new position onto the current primitive */ + /* Discard the hit if it is on a vertex/edge, and target a new position + * onto the current primitive */ } while((SXD_HIT_NONE(&hit) || HIT_ON_BOUNDARY(&hit)) && ++istep < nsteps); @@ -863,7 +863,7 @@ XD(scene_get_medium) cos_N_dir = fX(dot)(N, dir); /* Not too close and not roughly orthognonal */ - if(hit.distance > 1.e-6 || absf(cos_N_dir) > 1.e-1f) { + if(hit.distance > 1.e-6 && absf(cos_N_dir) > 1.e-1f) { const struct sdis_interface* interf; interf = scene_get_interface(scn, hit.prim.prim_id); medium = interface_get_medium diff --git a/src/test_sdis_solid_random_walk_robustness.c b/src/test_sdis_solid_random_walk_robustness.c @@ -0,0 +1,333 @@ +/* Copyright (C) 2016-2019 |Meso|Star> (contact@meso-star.com) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. */ + +#include "sdis.h" +#include "test_sdis_utils.h" + +#include <star/s3dut.h> +#include <rsys/math.h> + +#define Tfluid 350.0 /* Temperature of the fluid in Kelvin */ +#define Hcoef 1.0 /* Convection coefficient */ +#define Nreals 10000 /* #realisations */ +#define UNKNOWN -1 + +/******************************************************************************* + * Helper functions + ******************************************************************************/ +static void +compute_aabb + (const double* positions, + const size_t nvertices, + double lower[3], + double upper[3]) +{ + size_t i; + CHK(positions && nvertices && lower && upper); + + lower[0] = lower[1] = lower[2] = DBL_MAX; + upper[0] = upper[1] = upper[2] =-DBL_MAX; + + FOR_EACH(i, 0, nvertices) { + lower[0] = MMIN(lower[0], positions[i*3+0]); + lower[1] = MMIN(lower[1], positions[i*3+1]); + lower[2] = MMIN(lower[2], positions[i*3+2]); + upper[0] = MMAX(upper[0], positions[i*3+0]); + upper[1] = MMAX(upper[1], positions[i*3+1]); + upper[2] = MMAX(upper[2], positions[i*3+2]); + } +} + +static double +trilinear_temperature(const double pos[3]) +{ + const double a = 333; + const double b = 432; + const double c = 579; + double x, y, z; + CHK(pos[0] >= -10.0 && pos[0] <= 10.0); + CHK(pos[1] >= -10.0 && pos[1] <= 10.0); + CHK(pos[2] >= -10.0 && pos[2] <= 10.0); + + x = (pos[0] + 10.0) / 20.0; + y = (pos[1] + 10.0) / 20.0; + z = (pos[2] + 10.0) / 20.0; + + return a*x + b*y + c*z; +} + +/******************************************************************************* + * Geometry + ******************************************************************************/ +struct context { + struct s3dut_mesh_data msh; + struct sdis_interface* interf; +}; + +static void +get_indices(const size_t itri, size_t ids[3], void* context) +{ + const struct context* ctx = context; + CHK(ids && ctx && itri < ctx->msh.nprimitives); + ids[0] = ctx->msh.indices[itri*3+0]; + ids[2] = ctx->msh.indices[itri*3+1]; + ids[1] = ctx->msh.indices[itri*3+2]; +} + +static void +get_position(const size_t ivert, double pos[3], void* context) +{ + const struct context* ctx = context; + CHK(pos && ctx && ivert < ctx->msh.nvertices); + pos[0] = ctx->msh.positions[ivert*3+0]; + pos[1] = ctx->msh.positions[ivert*3+1]; + pos[2] = ctx->msh.positions[ivert*3+2]; +} + +static void +get_interface(const size_t itri, struct sdis_interface** bound, void* context) +{ + const struct context* ctx = context; + CHK(bound && ctx && itri < ctx->msh.nprimitives); + *bound = ctx->interf; +} + +/******************************************************************************* + * Interface + ******************************************************************************/ +struct interf { + double h; +}; + +static double +interface_get_temperature + (const struct sdis_interface_fragment* frag, struct sdis_data* data) +{ + const struct interf* interf; + CHK(data != NULL && frag != NULL); + interf = sdis_data_cget(data); + return interf->h > 0 ? UNKNOWN : trilinear_temperature(frag->P); +} + +static double +interface_get_convection_coef + (const struct sdis_interface_fragment* frag, struct sdis_data* data) +{ + CHK(frag != NULL); + return ((const struct interf*)sdis_data_cget(data))->h;; +} + +/******************************************************************************* + * Fluid + ******************************************************************************/ +static double +fluid_get_temperature + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + CHK(vtx != NULL); + (void)data; + return Tfluid; +} + +/******************************************************************************* + * Solid API + ******************************************************************************/ +struct solid { + double delta; + double cp; + double lambda; + double rho; + double temperature; +}; + +static double +solid_get_calorific_capacity + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + CHK(data != NULL && vtx != NULL); + return ((const struct solid*)sdis_data_cget(data))->cp; +} + +static double +solid_get_thermal_conductivity + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + CHK(data != NULL && vtx != NULL); + return ((const struct solid*)sdis_data_cget(data))->lambda; +} + +static double +solid_get_volumic_mass + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + CHK(data != NULL && vtx != NULL); + return ((const struct solid*)sdis_data_cget(data))->rho; +} + +static double +solid_get_delta + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + CHK(data != NULL && vtx != NULL); + return ((const struct solid*)sdis_data_cget(data))->delta; +} + +static double +solid_get_temperature + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + CHK(data != NULL && vtx != NULL); + return ((const struct solid*)sdis_data_cget(data))->temperature; +} + +/******************************************************************************* + * Main function + ******************************************************************************/ +int +main(int argc, char** argv) +{ + struct mem_allocator allocator; + struct s3dut_super_formula f0 = S3DUT_SUPER_FORMULA_NULL; + struct s3dut_super_formula f1 = S3DUT_SUPER_FORMULA_NULL; + struct s3dut_mesh* msh = NULL; + struct sdis_mc T = SDIS_MC_NULL; + struct sdis_device* dev = NULL; + struct sdis_data* data = NULL; + struct sdis_estimator* estimator = NULL; + struct sdis_medium* fluid = NULL; + struct sdis_medium* solid = NULL; + struct sdis_interface* interf = NULL; + struct sdis_scene* scn = NULL; + struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; + struct sdis_solid_shader solid_shader = SDIS_SOLID_SHADER_NULL; + struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; + struct interf* interf_param = NULL; + struct solid* solid_param = NULL; + struct context ctx; + double probe[3]; + double lower[3]; + double upper[3]; + double size[3]; + const double time[2] = {DBL_MAX, DBL_MAX}; + double ref; + size_t nreals; + size_t nfails; + (void)argc, (void)argv; + + OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator)); + OK(sdis_device_create(NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev)); + + /* Create the fluid medium */ + fluid_shader.temperature = fluid_get_temperature; + OK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid)); + + /* Setup the solid shader */ + solid_shader.calorific_capacity = solid_get_calorific_capacity; + solid_shader.thermal_conductivity = solid_get_thermal_conductivity; + solid_shader.volumic_mass = solid_get_volumic_mass; + solid_shader.delta_solid = solid_get_delta; + solid_shader.temperature = solid_get_temperature; + + /* Create the solid medium */ + OK(sdis_data_create + (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data)); + solid_param = sdis_data_get(data); + solid_param->delta = 0.01; + solid_param->lambda = 10; + solid_param->cp = 1.0; + solid_param->rho = 1.0; + solid_param->temperature = -1; + OK(sdis_solid_create(dev, &solid_shader, data, &solid)); + OK(sdis_data_ref_put(data)); + + /* Create the interface */ + OK(sdis_data_create + (dev, sizeof(struct interf), ALIGNOF(struct interf), NULL, &data)); + interf_param = sdis_data_get(data); + interf_param->h = 1; + interf_shader.convection_coef = interface_get_convection_coef; + interf_shader.front.temperature = interface_get_temperature; + OK(sdis_interface_create(dev, solid, fluid, &interf_shader, data, &interf)); + OK(sdis_data_ref_put(data)); + + /* Create the solid super shape */ + f0.A = 1; f0.B = 1; f0.M = 20; f0.N0 = 1; f0.N1 = 1; f0.N2 = 5; + f1.A = 1; f1.B = 1; f1.M = 7; f1.N0 = 1; f1.N1 = 2; f1.N2 = 5; + OK(s3dut_create_super_shape(&allocator, &f0, &f1, 1, 128, 64, &msh)); + OK(s3dut_mesh_get_data(msh, &ctx.msh)); + + compute_aabb(ctx.msh.positions, ctx.msh.nvertices, lower, upper); + + /* Create the scene */ + ctx.interf = interf; + OK(sdis_scene_create(dev, ctx.msh.nprimitives, get_indices, get_interface, + ctx.msh.nvertices, get_position, &ctx, &scn)); + /*dump_mesh(stdout, ctx.msh.positions, + ctx.msh.nvertices, ctx.msh.indices, ctx.msh.nprimitives);*/ + + /* Compute the delta of the solid random walk */ + size[0] = upper[0] - lower[0]; + size[1] = upper[1] - lower[1]; + size[2] = upper[2] - lower[2]; + solid_param->delta = MMIN(MMIN(size[0], size[1]), size[2]) / 20.0; + + /* Launch probe estimation */ + probe[0] = 0; + probe[1] = 0; + probe[2] = 0; + interf_param->h = 0; + OK(sdis_solve_probe(scn, Nreals, probe, time, 1.0, -1, 0, + SDIS_HEAT_PATH_FAILED, &estimator)); + + /* Print estimation results */ + OK(sdis_estimator_get_temperature(estimator, &T)); + OK(sdis_estimator_get_realisation_count(estimator, &nreals)); + OK(sdis_estimator_get_failure_count(estimator, &nfails)); + OK(sdis_estimator_ref_put(estimator)); + CHK(nfails < Nreals * 0.0001); + ref = trilinear_temperature(probe); + printf("T = %g ~ %g +/- %g; #failures = %lu / %lu\n", + ref, T.E, T.SE, (unsigned long)nfails, (unsigned long)Nreals); + CHK(eq_eps(T.E, ref, T.SE*3)); + + /* Launch medium integration */ + interf_param->h = 1; /* H delta T */ + OK(sdis_solve_medium(scn, Nreals, solid, time, 1.0, -1, 0, + SDIS_HEAT_PATH_FAILED, &estimator)); + /*dump_heat_paths(stdout, estimator);*/ + + /* Print estimation results */ + OK(sdis_estimator_get_temperature(estimator, &T)); + OK(sdis_estimator_get_realisation_count(estimator, &nreals)); + OK(sdis_estimator_get_failure_count(estimator, &nfails)); + OK(sdis_estimator_ref_put(estimator)); + printf("T ~ %g +/- %g; #failures = %lu / %lu\n", T.E, T.SE, + (unsigned long)nfails, (unsigned long)Nreals); + CHK(nfails < Nreals * 0.001); + + /* Release */ + OK(s3dut_mesh_ref_put(msh)); + OK(sdis_device_ref_put(dev)); + OK(sdis_medium_ref_put(fluid)); + OK(sdis_medium_ref_put(solid)); + OK(sdis_interface_ref_put(interf)); + OK(sdis_scene_ref_put(scn)); + + check_memory_allocator(&allocator); + mem_shutdown_proxy_allocator(&allocator); + CHK(mem_allocated_size() == 0); + return 0; +} + diff --git a/src/test_sdis_solve_probe3.c b/src/test_sdis_solve_probe3.c @@ -190,7 +190,7 @@ main(int argc, char** argv) double pos[3]; double time_range[2]; double ref; - const size_t N = 10000; + const size_t N = 100000; size_t ntris; size_t nverts; size_t nreals; @@ -312,7 +312,7 @@ main(int argc, char** argv) /* Check the results */ CHK(nfails + nreals == N); CHK(nfails < N/1000); - CHK(eq_eps(T.E, ref, 2*T.SE)); + CHK(eq_eps(T.E, ref, 3*T.SE)); /* Check green function */ OK(sdis_solve_probe_green_function(scn, N, pos, 1.0, -1, 0, &green)); diff --git a/src/test_sdis_volumic_power4_2d.c b/src/test_sdis_volumic_power4_2d.c @@ -22,7 +22,7 @@ #define Power 10000.0 #define H 50.0 #define LAMBDA 100.0 -#define DELTA (1.0/2.0) +#define DELTA 0.4/*(1.0/2.0)*/ #define N 10000 /*