diff options
Diffstat (limited to 'src/patch/patch_system.cc')
| -rw-r--r-- | src/patch/patch_system.cc | 246 |
1 files changed, 177 insertions, 69 deletions
diff --git a/src/patch/patch_system.cc b/src/patch/patch_system.cc index 575556f..c8ed7d0 100644 --- a/src/patch/patch_system.cc +++ b/src/patch/patch_system.cc @@ -23,10 +23,12 @@ // patch_system::patch_irho_isigma_of_gpn // patch_system::ghosted_patch_irho_isigma_of_gpn // +// patch_system::set_gridfn_to_constant // patch_system::gridfn_copy // patch_system::add_to_ghosted_gridfn // patch_system::gridfn_norms // patch_system::ghosted_gridfn_norms +// patch_system::integrate_gridfn // // patch_system::print_unknown_gridfn // patch_system::read_unknown_gridfn @@ -93,6 +95,9 @@ using jtutil::error_exit; // be 1; in general // N_overlap_points = 2*N_extend_points + 1 // delta_drho_dsigma = Grid spacing (both rho and sigma) in degrees. +// print_msg_flag = true to print CCTK_VInfo messages describing the +// construction process, +// false to skip this (and be silent except for errors) // patch_system::patch_system(fp origin_x_in, fp origin_y_in, fp origin_z_in, enum patch_system_type type_in, @@ -100,7 +105,8 @@ patch_system::patch_system(fp origin_x_in, fp origin_y_in, fp origin_z_in, fp delta_drho_dsigma, int min_gfn_in, int max_gfn_in, int ghosted_min_gfn_in, int ghosted_max_gfn_in, - int interp_handle, int interp_par_table_handle) + int interp_handle, int interp_par_table_handle, + bool print_msg_flag) : global_coords_(origin_x_in, origin_y_in, origin_z_in), type_(type_in), @@ -124,53 +130,63 @@ switch (type_in) { case full_sphere_patch_system: create_patches(patch_system_info::full_sphere::patch_info_array, - N_ghost_points, N_extend_points, - delta_drho_dsigma); + N_ghost_points, N_extend_points, delta_drho_dsigma, + print_msg_flag); setup_gridfn_storage(min_gfn_in, max_gfn_in, - ghosted_min_gfn_in, ghosted_max_gfn_in); + ghosted_min_gfn_in, ghosted_max_gfn_in, + print_msg_flag); interlink_full_sphere_patch_system(N_overlap_points, interp_handle, - interp_par_table_handle); + interp_par_table_handle, + print_msg_flag); break; case plus_z_hemisphere_patch_system: create_patches(patch_system_info::plus_z_hemisphere::patch_info_array, - N_ghost_points, N_extend_points, - delta_drho_dsigma); + N_ghost_points, N_extend_points, delta_drho_dsigma, + print_msg_flag); setup_gridfn_storage(min_gfn_in, max_gfn_in, - ghosted_min_gfn_in, ghosted_max_gfn_in); + ghosted_min_gfn_in, ghosted_max_gfn_in, + print_msg_flag); interlink_plus_z_hemisphere_patch_system(N_overlap_points, interp_handle, - interp_par_table_handle); + interp_par_table_handle, + print_msg_flag); break; case plus_xy_quadrant_patch_system: create_patches(patch_system_info::plus_xy_quadrant::patch_info_array, - N_ghost_points, N_extend_points, - delta_drho_dsigma); + N_ghost_points, N_extend_points, delta_drho_dsigma, + print_msg_flag); setup_gridfn_storage(min_gfn_in, max_gfn_in, - ghosted_min_gfn_in, ghosted_max_gfn_in); + ghosted_min_gfn_in, ghosted_max_gfn_in, + print_msg_flag); interlink_plus_xy_quadrant_patch_system(N_overlap_points, interp_handle, - interp_par_table_handle); + interp_par_table_handle, + print_msg_flag); break; case plus_xz_quadrant_patch_system: create_patches(patch_system_info::plus_xz_quadrant::patch_info_array, - N_ghost_points, N_extend_points, - delta_drho_dsigma); + N_ghost_points, N_extend_points, delta_drho_dsigma, + print_msg_flag); setup_gridfn_storage(min_gfn_in, max_gfn_in, - ghosted_min_gfn_in, ghosted_max_gfn_in); + ghosted_min_gfn_in, ghosted_max_gfn_in, + print_msg_flag); interlink_plus_xz_quadrant_patch_system(N_overlap_points, interp_handle, - interp_par_table_handle); + interp_par_table_handle, + print_msg_flag); break; case plus_xyz_octant_patch_system: create_patches(patch_system_info::plus_xyz_octant::patch_info_array, - N_ghost_points, N_extend_points, - delta_drho_dsigma); + N_ghost_points, N_extend_points, delta_drho_dsigma, + print_msg_flag); setup_gridfn_storage(min_gfn_in, max_gfn_in, - ghosted_min_gfn_in, ghosted_max_gfn_in); + ghosted_min_gfn_in, ghosted_max_gfn_in, + print_msg_flag); interlink_plus_xyz_octant_patch_system(N_overlap_points, interp_handle, - interp_par_table_handle); + interp_par_table_handle, + print_msg_flag); break; default: error_exit(ERROR_EXIT, @@ -213,23 +229,28 @@ delete gridfn_storage_; // void patch_system::create_patches(const struct patch_info patch_info_in[], int N_ghost_points, int N_extend_points, - fp delta_drho_dsigma) + fp delta_drho_dsigma, + bool print_msg_flag) { -CCTK_VInfo(CCTK_THORNSTRING, - " constructing %s patch system", - name_of_type(type())); -CCTK_VInfo(CCTK_THORNSTRING, - " at origin (%g,%g,%g)", - double(origin_x()), double(origin_y()), double(origin_z())); +if (print_msg_flag) + then { + CCTK_VInfo(CCTK_THORNSTRING, + " constructing %s patch system", + name_of_type(type())); + CCTK_VInfo(CCTK_THORNSTRING, + " at origin (%g,%g,%g)", + double(origin_x()), double(origin_y()), double(origin_z())); + } N_grid_points_ = 0; ghosted_N_grid_points_ = 0; for (int pn = 0 ; pn < N_patches() ; ++pn) { const struct patch_info& pi = patch_info_in[pn]; - CCTK_VInfo(CCTK_THORNSTRING, - " constructing %s patch", - pi.name); + if (print_msg_flag) + then CCTK_VInfo(CCTK_THORNSTRING, + " constructing %s patch", + pi.name); struct patch *p; switch (pi.ctype) @@ -322,7 +343,8 @@ ghosted_starting_gpn_[N_patches_] = ghosted_N_grid_points_; // void patch_system::setup_gridfn_storage (int min_gfn_in, int max_gfn_in, - int ghosted_min_gfn_in, int ghosted_max_gfn_in) + int ghosted_min_gfn_in, int ghosted_max_gfn_in, + bool print_msg_flag) { const int N_gridfns_in = jtutil::how_many_in_range(min_gfn_in, max_gfn_in); const int ghosted_N_gridfns_in @@ -334,14 +356,17 @@ const int ghosted_gfn_stride = ghosted_N_grid_points(); const int N_storage = gfn_stride * N_gridfns_in; const int ghosted_N_storage = ghosted_gfn_stride * ghosted_N_gridfns_in; -CCTK_VInfo(CCTK_THORNSTRING, " setting up gridfn storage"); -CCTK_VInfo(CCTK_THORNSTRING, - " gfn=[%d,%d] ghosted_gfn=[%d,%d]", - min_gfn_in, max_gfn_in, - ghosted_min_gfn_in, ghosted_max_gfn_in); -CCTK_VInfo(CCTK_THORNSTRING, - " N_grid_points()=%d ghosted_N_grid_points()=%d", - N_grid_points(), ghosted_N_grid_points()); +if (print_msg_flag) + then { + CCTK_VInfo(CCTK_THORNSTRING, " setting up gridfn storage"); + CCTK_VInfo(CCTK_THORNSTRING, + " gfn=[%d,%d] ghosted_gfn=[%d,%d]", + min_gfn_in, max_gfn_in, + ghosted_min_gfn_in, ghosted_max_gfn_in); + CCTK_VInfo(CCTK_THORNSTRING, + " N_grid_points()=%d ghosted_N_grid_points()=%d", + N_grid_points(), ghosted_N_grid_points()); + } // storage arrays for all gridfns gridfn_storage_ = new fp[N_storage]; @@ -370,8 +395,11 @@ ghosted_gridfn_storage_ = new fp[ghosted_N_storage]; p.setup_gridfn_storage(gridfn_pars, ghosted_gridfn_pars); } -CCTK_VInfo(CCTK_THORNSTRING, - " checking that storage is partitioned properly"); +if (print_msg_flag) + then { + CCTK_VInfo(CCTK_THORNSTRING, + " checking that storage is partitioned properly"); + } // check to make sure storage for distinct gridfns // forms a partition of the overall storage array @@ -508,9 +536,12 @@ const patch& plast = ith_patch(N_patches()-1); // void patch_system::interlink_full_sphere_patch_system (int N_overlap_points, - int interp_handle, int interp_par_table_handle) + int interp_handle, int interp_par_table_handle, + bool print_msg_flag) { -CCTK_VInfo(CCTK_THORNSTRING, " interlinking full sphere patch system"); +if (print_msg_flag) + then CCTK_VInfo(CCTK_THORNSTRING, + " interlinking full sphere patch system"); patch& pz = ith_patch(patch_number_of_name("+z")); patch& px = ith_patch(patch_number_of_name("+x")); @@ -520,7 +551,8 @@ patch& my = ith_patch(patch_number_of_name("-y")); patch& mz = ith_patch(patch_number_of_name("-z")); // create the ghost zones -CCTK_VInfo(CCTK_THORNSTRING, " creating ghost zones"); +if (print_msg_flag) + then CCTK_VInfo(CCTK_THORNSTRING, " creating ghost zones"); create_interpatch_ghost_zones(pz, px, N_overlap_points); create_interpatch_ghost_zones(pz, py, N_overlap_points); create_interpatch_ghost_zones(pz, mx, N_overlap_points); @@ -535,7 +567,8 @@ create_interpatch_ghost_zones(mz, mx, N_overlap_points); create_interpatch_ghost_zones(mz, my, N_overlap_points); // finish setting up the interpatch ghost zones -CCTK_VInfo(CCTK_THORNSTRING, " finishing interpatch setup"); +if (print_msg_flag) + then CCTK_VInfo(CCTK_THORNSTRING, " finishing interpatch setup"); finish_interpatch_setup(pz, px, N_overlap_points, interp_handle, interp_par_table_handle); @@ -584,9 +617,12 @@ assert_all_ghost_zones_fully_setup(); // void patch_system::interlink_plus_z_hemisphere_patch_system (int N_overlap_points, - int interp_handle, int interp_par_table_handle) + int interp_handle, int interp_par_table_handle, + bool print_msg_flag) { -CCTK_VInfo(CCTK_THORNSTRING, " interlinking +z hemisphere patch system"); +if (print_msg_flag) + then CCTK_VInfo(CCTK_THORNSTRING, + " interlinking +z hemisphere patch system"); patch& pz = ith_patch(patch_number_of_name("+z")); patch& px = ith_patch(patch_number_of_name("+x")); @@ -595,7 +631,8 @@ patch& mx = ith_patch(patch_number_of_name("-x")); patch& my = ith_patch(patch_number_of_name("-y")); // create the ghost zones -CCTK_VInfo(CCTK_THORNSTRING, " creating ghost zones"); +if (print_msg_flag) + then CCTK_VInfo(CCTK_THORNSTRING, " creating ghost zones"); create_interpatch_ghost_zones(pz, px, N_overlap_points); create_interpatch_ghost_zones(pz, py, N_overlap_points); create_interpatch_ghost_zones(pz, mx, N_overlap_points); @@ -610,7 +647,8 @@ mx.create_mirror_symmetry_ghost_zone(mx.min_rho_patch_edge()); my.create_mirror_symmetry_ghost_zone(my.min_rho_patch_edge()); // finish setting up the interpatch ghost zones -CCTK_VInfo(CCTK_THORNSTRING, " finishing interpatch setup"); +if (print_msg_flag) + then CCTK_VInfo(CCTK_THORNSTRING, " finishing interpatch setup"); finish_interpatch_setup(pz, px, N_overlap_points, interp_handle, interp_par_table_handle); @@ -647,9 +685,12 @@ assert_all_ghost_zones_fully_setup(); // void patch_system::interlink_plus_xy_quadrant_patch_system (int N_overlap_points, - int interp_handle, int interp_par_table_handle) + int interp_handle, int interp_par_table_handle, + bool print_msg_flag) { -CCTK_VInfo(CCTK_THORNSTRING, " interlinking +xy quadrant patch system"); +if (print_msg_flag) + then CCTK_VInfo(CCTK_THORNSTRING, + " interlinking +xy quadrant patch system"); patch& pz = ith_patch(patch_number_of_name("+z")); patch& px = ith_patch(patch_number_of_name("+x")); @@ -657,7 +698,8 @@ patch& py = ith_patch(patch_number_of_name("+y")); patch& mz = ith_patch(patch_number_of_name("-z")); // create the ghost zones -CCTK_VInfo(CCTK_THORNSTRING, " creating ghost zones"); +if (print_msg_flag) + then CCTK_VInfo(CCTK_THORNSTRING, " creating ghost zones"); create_interpatch_ghost_zones(pz, px, N_overlap_points); create_interpatch_ghost_zones(pz, py, N_overlap_points); create_interpatch_ghost_zones(px, py, N_overlap_points); @@ -674,7 +716,8 @@ create_periodic_symmetry_ghost_zones(mz.max_rho_patch_edge(), true); // finish setting up the interpatch ghost zones -CCTK_VInfo(CCTK_THORNSTRING, " finishing interpatch setup"); +if (print_msg_flag) + then CCTK_VInfo(CCTK_THORNSTRING, " finishing interpatch setup"); finish_interpatch_setup(pz, px, N_overlap_points, interp_handle, interp_par_table_handle); @@ -702,9 +745,12 @@ assert_all_ghost_zones_fully_setup(); // void patch_system::interlink_plus_xz_quadrant_patch_system (int N_overlap_points, - int interp_handle, int interp_par_table_handle) + int interp_handle, int interp_par_table_handle, + bool print_msg_flag) { -CCTK_VInfo(CCTK_THORNSTRING, " interlinking +xz quadrant patch system"); +if (print_msg_flag) + then CCTK_VInfo(CCTK_THORNSTRING, + " interlinking +xz quadrant patch system"); patch& pz = ith_patch(patch_number_of_name("+z")); patch& px = ith_patch(patch_number_of_name("+x")); @@ -712,7 +758,8 @@ patch& py = ith_patch(patch_number_of_name("+y")); patch& my = ith_patch(patch_number_of_name("-y")); // create the ghost zones -CCTK_VInfo(CCTK_THORNSTRING, " creating ghost zones"); +if (print_msg_flag) + then CCTK_VInfo(CCTK_THORNSTRING, " creating ghost zones"); create_interpatch_ghost_zones(pz, px, N_overlap_points); create_interpatch_ghost_zones(pz, py, N_overlap_points); create_interpatch_ghost_zones(pz, my, N_overlap_points); @@ -729,7 +776,8 @@ create_periodic_symmetry_ghost_zones(py.max_sigma_patch_edge(), false); // finish setting up the interpatch ghost zones -CCTK_VInfo(CCTK_THORNSTRING, " finishing interpatch setup"); +if (print_msg_flag) + then CCTK_VInfo(CCTK_THORNSTRING, " finishing interpatch setup"); finish_interpatch_setup(pz, px, N_overlap_points, interp_handle, interp_par_table_handle); @@ -757,16 +805,20 @@ assert_all_ghost_zones_fully_setup(); // void patch_system::interlink_plus_xyz_octant_patch_system (int N_overlap_points, - int interp_handle, int interp_par_table_handle) + int interp_handle, int interp_par_table_handle, + bool print_msg_flag) { -CCTK_VInfo(CCTK_THORNSTRING, " interlinking +xyz octant patch system"); +if (print_msg_flag) + then CCTK_VInfo(CCTK_THORNSTRING, + " interlinking +xyz octant patch system"); patch& pz = ith_patch(patch_number_of_name("+z")); patch& px = ith_patch(patch_number_of_name("+x")); patch& py = ith_patch(patch_number_of_name("+y")); // create the ghost zones -CCTK_VInfo(CCTK_THORNSTRING, " creating ghost zones"); +if (print_msg_flag) + then CCTK_VInfo(CCTK_THORNSTRING, " creating ghost zones"); create_interpatch_ghost_zones(pz, px, N_overlap_points); create_interpatch_ghost_zones(pz, py, N_overlap_points); create_interpatch_ghost_zones(px, py, N_overlap_points); @@ -780,7 +832,8 @@ create_periodic_symmetry_ghost_zones(px.min_sigma_patch_edge(), true); // finish setting up the interpatch ghost zones -CCTK_VInfo(CCTK_THORNSTRING, " finishing interpatch setup"); +if (print_msg_flag) + then CCTK_VInfo(CCTK_THORNSTRING, " finishing interpatch setup"); finish_interpatch_setup(pz, px, N_overlap_points, interp_handle, interp_par_table_handle); @@ -990,7 +1043,8 @@ error_exit(ERROR_EXIT, // design would be to return this via a patch*& argument, but the design // used here seems slightly cleaner to use in practice.) // -patch& patch_system::patch_irho_isigma_of_gpn(int gpn, int& irho, int& isigma) +const patch& + patch_system::patch_irho_isigma_of_gpn(int gpn, int& irho, int& isigma) const { assert( gpn >= 0 ); @@ -1001,7 +1055,7 @@ assert( gpn < N_grid_points() ); // n.b. [pn+1] is ok since starting_gpn_[] has size N_patches()+1 if ((gpn >= starting_gpn_[pn]) && (gpn < starting_gpn_[pn+1])) then { - patch& p = ith_patch(pn); + const patch& p = ith_patch(pn); const int gpn_in_patch = gpn - starting_gpn_[pn]; assert( gpn_in_patch >= 0 ); assert( gpn_in_patch < p.N_grid_points() ); @@ -1034,7 +1088,8 @@ error_exit(PANIC_EXIT, // design would be to return this via a patch*& argument, but the design // used here seems slightly cleaner to use in practice.) // -patch& patch_system::ghosted_patch_irho_isigma_of_gpn +const patch& + patch_system::ghosted_patch_irho_isigma_of_gpn (int gpn, int& irho, int& isigma) const { @@ -1048,7 +1103,7 @@ assert( gpn < ghosted_N_grid_points() ); if ( (gpn >= ghosted_starting_gpn_[pn]) && (gpn < ghosted_starting_gpn_[pn+1])) then { - patch& p = ith_patch(pn); + const patch& p = ith_patch(pn); const int gpn_in_patch = gpn - ghosted_starting_gpn_[pn]; assert( gpn_in_patch >= 0 ); assert( gpn_in_patch < p.ghosted_N_grid_points() ); @@ -1071,6 +1126,28 @@ error_exit(PANIC_EXIT, //****************************************************************************** // +// This function sets a (nominal-grid) gridfn to a constant value. +// +void patch_system::set_gridfn_to_constant(fp a, int dst_gfn) +{ + for (int pn = 0 ; pn < N_patches() ; ++pn) + { + patch& p = ith_patch(pn); + for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho) + { + for (int isigma = p.min_isigma() ; + isigma <= p.max_isigma() ; + ++isigma) + { + p.gridfn(dst_gfn, irho,isigma) = a; + } + } + } +} + +//****************************************************************************** + +// // This function copies one (nominal-grid) gridfn to another. // void patch_system::gridfn_copy(int src_gfn, int dst_gfn) @@ -1120,6 +1197,7 @@ void patch_system::add_to_ghosted_gridfn(fp delta, int ghosted_dst_gfn) // This function computes norms of a nominal-grid gridfn. // void patch_system::gridfn_norms(int src_gfn, jtutil::norm<fp>& norms) + const { if (! is_valid_gfn(src_gfn)) then error_exit(ERROR_EXIT, @@ -1130,7 +1208,7 @@ norms.reset(); for (int pn = 0 ; pn < N_patches() ; ++pn) { - patch& p = ith_patch(pn); + const patch& p = ith_patch(pn); for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho) { for (int isigma = p.min_isigma() ; @@ -1150,6 +1228,7 @@ norms.reset(); // void patch_system::ghosted_gridfn_norms(int ghosted_src_gfn, jtutil::norm<fp>& norms) + const { if (! is_valid_ghosted_gfn(ghosted_src_gfn)) then error_exit(ERROR_EXIT, @@ -1159,7 +1238,7 @@ norms.reset(); for (int pn = 0 ; pn < N_patches() ; ++pn) { - patch& p = ith_patch(pn); + const patch& p = ith_patch(pn); for (int irho = p.ghosted_min_irho() ; irho <= p.ghosted_max_irho() ; ++irho) @@ -1175,6 +1254,35 @@ norms.reset(); } //****************************************************************************** + +// +// This function computes an approximation to the (surface) integral of +// a gridfn over the whole patch system, +// area_weighting_flag ? $\int f(\rho,\sigma) \, dA$ +// : $\int f(\rho,\sigma) \, d\rho \, d\sigma$ +// where $f(\rho,\sigma)$ is the gridfn and $dA$ is the area element in +// $(\rho,sigma)$ coordinates. +// +// The integration scheme is selected based on the method argument. +// +// Bugs: +// The way the integration coefficients are computed is somewhat inefficient. +// +fp patch_system::integrate_gridfn(int src_gfn, + bool area_weighting_flag, + enum patch::integration_method method) + const +{ +fp integral = 0.0; + for (int pn = 0 ; pn < N_patches() ; ++pn) + { + const patch& p = ith_patch(pn); + integral += p.integrate_gridfn(src_gfn, area_weighting_flag, method); + } +return integral; +} + +//****************************************************************************** //****************************************************************************** //****************************************************************************** @@ -1483,7 +1591,7 @@ global_min_ym_ = +INT_MAX; global_max_ym_ = -INT_MAX; for (int pn = 0 ; pn < N_patches() ; ++pn) { - patch& p = ith_patch(pn); + const patch& p = ith_patch(pn); // n.b. these loops must use _int_ variables for the loop // to terminate! for (int want_min = false ; want_min <= true ; ++want_min) |