diff --git a/docs/documentation/case.md b/docs/documentation/case.md index b92e8f945e..92661ef182 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -312,6 +312,8 @@ This is enabled by adding ``'elliptic_smoothing': "T",`` and ``'elliptic_smoothi | Parameter | Type | Description | | ---: | :----: | :--- | | `num_ibs` | Integer | Number of immersed boundary patches | +| `num_particle_beds` | Integer | Number of particle bed specifications to generate immersed boundary patches from | +| `ib_neighborhood_radius` | Integer | Parameter that controls the neighborhood size for IB detection. | | `geometry` | Integer | Geometry configuration of the patch.| | `x[y,z]_centroid` | Real | Centroid of the applied geometry in the [x,y,z]-direction. | | `length_x[y,z]` | Real | Length, if applicable, in the [x,y,z]-direction. | @@ -373,6 +375,8 @@ Additional details on this specification can be found in [NACA airfoil](https:// - `ib_coefficient_of_friction` is the coefficient of friction used in IB collisions. +- `ib_neighborhood_radius` controls the size of the neighborhood size. This value defaults to 1, which indicates that any given rank is aware of IB's up to 1 ranks away. This parameter is required to strong-scale a case when IB's eventually grow to be larger than one full processor domain wide. + ### 5. Fluid Material's {#sec-fluid-materials} | Parameter | Type | Description | @@ -643,7 +647,7 @@ To restart the simulation from $k$-th time step, see @ref running "Restarting Ca | `alpha_wrt(i)` | Logical | Add the volume fraction of fluid $i$ to the database | | `gamma_wrt` | Logical | Add the specific heat ratio function to the database | | `heat_ratio_wrt` | Logical | Add the specific heat ratio to the database | -| `ib_state_wrt` | Logical | Write IB state and loads to a datafile at each time step | +| `ib_state_wrt` | Logical | Parameter to handle writing IB state on saves and outputting the state as a point mesh to SILO files. | | `pi_inf_wrt` | Logical | Add the liquid stiffness function to the database | | `pres_inf_wrt` | Logical | Add the liquid stiffness to the formatted database | | `c_wrt` | Logical | Add the sound speed to the database | @@ -711,7 +715,7 @@ If `file_per_process` is true, then pre_process, simulation, and post_process mu - `probe_wrt` activates the output of state variables at coordinates specified by `probe(i)%[x;y,z]`. -- `ib_state_wrt` activates the output of data specified by patch_ib(i)%force(:) (and torque, vel, angular_vel, angles, [x,y,z]_centroid) into a single binary datafile for all IBs at all timesteps. During post_processing, this file is converted into separate time histories for each IB. +- `ib_state_wrt` is used to trigger post-processing of the IB state to be written out as a point mesh in the SILO files. When no IBs are moving, it also triggers force and torque calculation so that those values may be written to the output state files. - `output_partial_domain` activates the output of part of the domain specified by `[x,y,z]_output%%beg` and `[x,y,z]_output%%end`. This is useful for large domains where only a portion of the domain is of interest. diff --git a/docs/module_categories.json b/docs/module_categories.json index ef4a7d726e..21de015c71 100644 --- a/docs/module_categories.json +++ b/docs/module_categories.json @@ -38,6 +38,7 @@ "m_compute_cbc", "m_boundary_common", "m_ibm", + "m_particle_bed", "m_igr", "m_ib_patches", "m_compute_levelset" diff --git a/examples/2D_mibm_particle_bed/case.py b/examples/2D_mibm_particle_bed/case.py new file mode 100644 index 0000000000..5db1e707c7 --- /dev/null +++ b/examples/2D_mibm_particle_bed/case.py @@ -0,0 +1,136 @@ +import json +import math + +# 2D shock wave interacting with a bed of 20 free-floating circular particles. + +gam_a = 1.4 + +# Shock parameters (Mach 1.5) +mach_number = 1.5 +pre_shock_pressure = 1 +pre_shock_density = 1.4 +pre_shock_speed = 0.0 +post_shock_pressure = 2.4583 +post_shock_density = 2.6069 +post_shock_speed = 0.6944 + +domain_size = 4.0 +wave_front = -1.5 + +total_time = 1.5 +num_time_steps = 2000 +dt = float(total_time / num_time_steps) +num_saves = 100 +steps_to_save = int(num_time_steps / num_saves) + +# Soft-sphere collision parameters (from 3D_mibm_sphere_head_on_collision) +collision_time = 20.0 * dt + +# Particle bed parameters +bed_x = 0.5 +bed_y = 0.0 +bed_lx = 2.0 +bed_ly = 3.5 +particle_radius = 0.15 +particle_mass = 0.25 +particle_min_spacing = 0.05 + +print( + json.dumps( + { + # Logistics + "run_time_info": "T", + # Computational Domain Parameters + "x_domain%beg": -domain_size * 0.5, + "x_domain%end": domain_size * 0.5, + "y_domain%beg": -domain_size * 0.5, + "y_domain%end": domain_size * 0.5, + "cyl_coord": "F", + "m": 256, + "n": 256, + "p": 0, + "dt": dt, + "t_step_start": 0, + "t_step_stop": num_time_steps, + "t_step_save": steps_to_save, + # Simulation Algorithm Parameters + "num_patches": 2, + "model_eqns": 2, + "alt_soundspeed": "F", + "num_fluids": 1, + "mpp_lim": "F", + "mixture_err": "T", + "time_stepper": 3, + "weno_order": 5, + "weno_eps": 1.0e-16, + "weno_Re_flux": "T", + "weno_avg": "T", + "avg_state": 2, + "mapped_weno": "T", + "null_weights": "F", + "mp_weno": "T", + "riemann_solver": 2, + "wave_speeds": 1, + "bc_x%beg": -17, + "bc_x%end": -8, + "bc_y%beg": -15, + "bc_y%end": -15, + # Immersed boundaries — all circles come from the particle bed + "ib": "T", + "num_ibs": 0, + "viscous": "T", + # Collision model (soft-sphere, from 3D_mibm_sphere_head_on_collision) + "collision_model": 1, + "coefficient_of_restitution": 0.9, + "collision_time": collision_time, + "ib_coefficient_of_friction": 0.1, + # Particle bed: 20 free-floating circles placed randomly in region + "num_particle_beds": 1, + "particle_bed(1)%x_centroid": bed_x, + "particle_bed(1)%y_centroid": bed_y, + "particle_bed(1)%z_centroid": 0.0, + "particle_bed(1)%length_x": bed_lx, + "particle_bed(1)%length_y": bed_ly, + "particle_bed(1)%length_z": 0.0, + "particle_bed(1)%num_particles": 20, + "particle_bed(1)%radius": particle_radius, + "particle_bed(1)%mass": particle_mass, + "particle_bed(1)%min_spacing": particle_min_spacing, + "particle_bed(1)%moving_ibm": 2, + "particle_bed(1)%seed": 42, + # Output + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "E_wrt": "T", + "ib_state_wrt": "F", + "parallel_io": "T", + # IC Patch 1: post-shock region (left of wave front) + "patch_icpp(1)%geometry": 3, + "patch_icpp(1)%x_centroid": 0.5 * wave_front - 0.25 * domain_size, + "patch_icpp(1)%y_centroid": 0.0, + "patch_icpp(1)%length_x": 0.5 * domain_size + wave_front, + "patch_icpp(1)%length_y": domain_size, + "patch_icpp(1)%vel(1)": post_shock_speed, + "patch_icpp(1)%vel(2)": 0.0, + "patch_icpp(1)%pres": post_shock_pressure, + "patch_icpp(1)%alpha_rho(1)": post_shock_density, + "patch_icpp(1)%alpha(1)": 1.0, + # IC Patch 2: pre-shock region (right of wave front) + "patch_icpp(2)%geometry": 3, + "patch_icpp(2)%x_centroid": 0.5 * wave_front + 0.25 * domain_size, + "patch_icpp(2)%y_centroid": 0.0, + "patch_icpp(2)%length_x": 0.5 * domain_size - wave_front, + "patch_icpp(2)%length_y": domain_size, + "patch_icpp(2)%vel(1)": pre_shock_speed, + "patch_icpp(2)%vel(2)": 0.0, + "patch_icpp(2)%pres": pre_shock_pressure, + "patch_icpp(2)%alpha_rho(1)": pre_shock_density, + "patch_icpp(2)%alpha(1)": 1.0, + # Fluid properties: air + "fluid_pp(1)%gamma": 1.0 / (gam_a - 1.0), + "fluid_pp(1)%pi_inf": 0, + "fluid_pp(1)%Re(1)": 2500000, + } + ) +) diff --git a/examples/2D_mibm_shock_cylinder/case.py b/examples/2D_mibm_shock_cylinder/case.py index 69aca0cc9b..9fbc86829d 100644 --- a/examples/2D_mibm_shock_cylinder/case.py +++ b/examples/2D_mibm_shock_cylinder/case.py @@ -87,7 +87,7 @@ "precision": 2, "prim_vars_wrt": "T", "E_wrt": "T", - "ib_state_wrt": "T", + "ib_state_wrt": "F", "parallel_io": "T", # Patch: Constant Tube filled with air # Specify the cylindrical air tube grid geometry diff --git a/src/common/include/case.fpp b/src/common/include/case.fpp index aa0e0637b9..8f5fc4777b 100644 --- a/src/common/include/case.fpp +++ b/src/common/include/case.fpp @@ -4,10 +4,8 @@ ! For pre-process. #:def analytical() - #:enddef ! For moving immersed boundaries in simulation #:def mib_analytical() - #:enddef diff --git a/src/common/m_constants.fpp b/src/common/m_constants.fpp index 77e369c7f5..34332850fa 100644 --- a/src/common/m_constants.fpp +++ b/src/common/m_constants.fpp @@ -16,15 +16,19 @@ module m_constants real(wp), parameter :: verysmall = 1.e-12_wp !< Very small number !> Radius cutoff to avoid division by zero for 3D spherical harmonic patch (geometry 14) real(wp), parameter :: small_radius = 1.e-32_wp - integer, parameter :: num_stcls_min = 5 !< Minimum # of stencils - integer, parameter :: path_len = 400 !< Maximum path length - integer, parameter :: name_len = 50 !< Maximum name length - integer, parameter :: dflt_int = -100 !< Default integer value - integer, parameter :: fourier_rings = 5 !< Fourier filter ring limit - integer, parameter :: num_fluids_max = 10 !< Maximum number of fluids in the simulation - integer, parameter :: num_probes_max = 10 !< Maximum number of flow probes in the simulation - integer, parameter :: num_patches_max = 10 !< Maximum number of IC patches - integer, parameter :: num_ib_patches_max = 50000 !< Maximum number of immersed boundary patches (patch_ib) + integer, parameter :: num_stcls_min = 5 !< Minimum # of stencils + integer, parameter :: path_len = 400 !< Maximum path length + integer, parameter :: name_len = 50 !< Maximum name length + integer, parameter :: dflt_int = -100 !< Default integer value + integer, parameter :: fourier_rings = 5 !< Fourier filter ring limit + integer, parameter :: num_fluids_max = 10 !< Maximum number of fluids in the simulation + integer, parameter :: num_probes_max = 10 !< Maximum number of flow probes in the simulation + integer, parameter :: num_patches_max = 10 !< Maximum number of IC patches + integer, parameter :: num_ib_patches_max = 2050000 !< Maximum number of immersed boundary patches (patch_ib) + !> Max patches readable from the namelist; patch_ib grows to num_ib_patches_max only when particle beds are used + integer, parameter :: num_ib_patches_max_namelist = 50000 + integer, parameter :: num_local_ibs_max = 2000 !< Maximum number of immersed boundary patches (patch_ib) + integer, parameter :: num_particle_beds_max = 10 !< Maximum number of particle bed patch specifications integer, parameter :: num_bc_patches_max = 10 !< Maximum number of boundary condition patches integer, parameter :: max_2d_fourier_modes = 10 !< Max Fourier mode index for 2D modal patch (geometry 13) integer, parameter :: max_sph_harm_degree = 5 !< Max degree L for 3D spherical harmonic patch (geometry 14) diff --git a/src/common/m_derived_types.fpp b/src/common/m_derived_types.fpp index 6ee2e836a5..c4e5873ead 100644 --- a/src/common/m_derived_types.fpp +++ b/src/common/m_derived_types.fpp @@ -202,12 +202,10 @@ module m_derived_types type :: t_model_array ! Original CPU-side fields (unchanged) - type(t_model), allocatable :: model !< STL/OBJ geometry model - real(wp), allocatable, dimension(:,:,:) :: boundary_v !< Boundary vertices - real(wp), allocatable, dimension(:,:) :: interpolated_boundary_v !< Interpolated boundary vertices - integer :: boundary_edge_count !< Number of boundary edges - integer :: total_vertices !< Total vertex count - integer :: interpolate !< Interpolation flag + type(t_model), allocatable :: model !< STL/OBJ geometry model + real(wp), allocatable, dimension(:,:,:) :: boundary_v !< Boundary vertices + integer :: boundary_edge_count !< Number of boundary edges + integer :: total_vertices !< Total vertex count ! GPU-friendly flattened arrays integer :: ntrs !< Copy of model%ntrs @@ -272,9 +270,10 @@ module m_derived_types end type ic_patch_parameters type ib_patch_parameters - integer :: geometry !< Type of geometry for the patch + integer :: gbl_patch_id real(wp) :: x_centroid, y_centroid, z_centroid !< Geometric center coordinates of the patch + !> Centroid locations of intermediate steps in the time_stepper module real(wp) :: step_x_centroid, step_y_centroid, step_z_centroid real(wp), dimension(1:3) :: centroid_offset !< offset of center of mass from computed cell center for odd-shaped IBs diff --git a/src/common/m_model.fpp b/src/common/m_model.fpp index 0dab036f9b..a526eb1266 100644 --- a/src/common/m_model.fpp +++ b/src/common/m_model.fpp @@ -983,12 +983,13 @@ contains dx_local = minval(dx); dy_local = minval(dy) if (p /= 0) dz_local = minval(dz) + num_gbl_ibs = num_ibs allocate (stl_bounding_boxes(num_ibs,1:3,1:3)) do patch_id = 1, num_ibs if (patch_ib(patch_id)%geometry == 5 .or. patch_ib(patch_id)%geometry == 12) then allocate (models(patch_id)%model) - print *, " * Reading model: " // trim(patch_ib(patch_id)%model_filepath) + if (proc_rank == 0) print *, " * Reading model: " // trim(patch_ib(patch_id)%model_filepath) model = f_model_read(patch_ib(patch_id)%model_filepath) params%scale(:) = patch_ib(patch_id)%model_scale(:) @@ -1001,9 +1002,7 @@ contains params%scale(:) = 1._wp end if - if (proc_rank == 0) then - print *, " * Transforming model." - end if + if (proc_rank == 0) print *, " * Transforming model." ! Get the model center before transforming the model bbox_old = f_create_bbox(model) diff --git a/src/common/m_mpi_common.fpp b/src/common/m_mpi_common.fpp index 7d4b92705c..9cfdd82e14 100644 --- a/src/common/m_mpi_common.fpp +++ b/src/common/m_mpi_common.fpp @@ -94,6 +94,8 @@ contains proc_rank = 0 #endif + $:GPU_UPDATE(device='[num_procs, proc_rank]') + end subroutine s_mpi_initialize !> Set up MPI I/O data views and variable pointers for parallel file output. diff --git a/src/post_process/m_data_output.fpp b/src/post_process/m_data_output.fpp index aece31ff8d..3a3a7bb10c 100644 --- a/src/post_process/m_data_output.fpp +++ b/src/post_process/m_data_output.fpp @@ -1396,48 +1396,42 @@ contains end do close (file_unit) - end if - - call MPI_BCAST(ib_data, nBodies*NFIELDS_PER_IB, mpi_p, 0, MPI_COMM_WORLD, ierr) - - do i = 1, nBodies - force_x(i) = ib_data(i, 2); force_y(i) = ib_data(i, 3); force_z(i) = ib_data(i, 4) - torque_x(i) = ib_data(i, 5); torque_y(i) = ib_data(i, 6); torque_z(i) = ib_data(i, 7) - vel_x(i) = ib_data(i, 8); vel_y(i) = ib_data(i, 9); vel_z(i) = ib_data(i, 10) - omega_x(i) = ib_data(i, 11); omega_y(i) = ib_data(i, 12); omega_z(i) = ib_data(i, 13) - angle_x(i) = ib_data(i, 14); angle_y(i) = ib_data(i, 15); angle_z(i) = ib_data(i, 16) - px(i) = ib_data(i, 17); py(i) = ib_data(i, 18); pz(i) = ib_data(i, 19) - ib_diameter(i) = ib_data(i, 20)*2.0_wp - end do - if (proc_rank == 0) then - do i = 1, num_procs - write (meshnames(i), '(A,I0,A,I0,A)') '../p', i - 1, '/', t_step, '.silo:ib_bodies' - meshtypes(i) = DB_POINTMESH + do i = 1, nBodies + force_x(i) = ib_data(i, 2); force_y(i) = ib_data(i, 3); force_z(i) = ib_data(i, 4) + torque_x(i) = ib_data(i, 5); torque_y(i) = ib_data(i, 6); torque_z(i) = ib_data(i, 7) + vel_x(i) = ib_data(i, 8); vel_y(i) = ib_data(i, 9); vel_z(i) = ib_data(i, 10) + omega_x(i) = ib_data(i, 11); omega_y(i) = ib_data(i, 12); omega_z(i) = ib_data(i, 13) + angle_x(i) = ib_data(i, 14); angle_y(i) = ib_data(i, 15); angle_z(i) = ib_data(i, 16) + px(i) = ib_data(i, 17); py(i) = ib_data(i, 18); pz(i) = ib_data(i, 19) + ib_diameter(i) = ib_data(i, 20)*2.0_wp end do + + write (meshnames(1), '(A,I0,A)') '../p0/', t_step, '.silo:ib_bodies' + meshtypes(1) = DB_POINTMESH err = DBSET2DSTRLEN(len(meshnames(1))) - err = DBPUTMMESH(dbroot, 'ib_bodies', 16, num_procs, meshnames, len_trim(meshnames), meshtypes, DB_F77NULL, ierr) + err = DBPUTMMESH(dbroot, 'ib_bodies', 16, 1, meshnames, len_trim(meshnames), meshtypes, DB_F77NULL, ierr) + + err = DBPUTPM(dbfile, 'ib_bodies', 9, 3, px, py, pz, nBodies, DB_DOUBLE, DB_F77NULL, ierr) + + call s_write_ib_variable('ib_force_x', t_step, force_x, nBodies) + call s_write_ib_variable('ib_force_y', t_step, force_y, nBodies) + call s_write_ib_variable('ib_force_z', t_step, force_z, nBodies) + call s_write_ib_variable('ib_torque_x', t_step, torque_x, nBodies) + call s_write_ib_variable('ib_torque_y', t_step, torque_y, nBodies) + call s_write_ib_variable('ib_torque_z', t_step, torque_z, nBodies) + call s_write_ib_variable('ib_vel_x', t_step, vel_x, nBodies) + call s_write_ib_variable('ib_vel_y', t_step, vel_y, nBodies) + call s_write_ib_variable('ib_vel_z', t_step, vel_z, nBodies) + call s_write_ib_variable('ib_omega_x', t_step, omega_x, nBodies) + call s_write_ib_variable('ib_omega_y', t_step, omega_y, nBodies) + call s_write_ib_variable('ib_omega_z', t_step, omega_z, nBodies) + call s_write_ib_variable('ib_angle_x', t_step, angle_x, nBodies) + call s_write_ib_variable('ib_angle_y', t_step, angle_y, nBodies) + call s_write_ib_variable('ib_angle_z', t_step, angle_z, nBodies) + call s_write_ib_variable('ib_diameter', t_step, ib_diameter, nBodies) end if - err = DBPUTPM(dbfile, 'ib_bodies', 9, 3, px, py, pz, nBodies, DB_DOUBLE, DB_F77NULL, ierr) - - call s_write_ib_variable('ib_force_x', t_step, force_x, nBodies) - call s_write_ib_variable('ib_force_y', t_step, force_y, nBodies) - call s_write_ib_variable('ib_force_z', t_step, force_z, nBodies) - call s_write_ib_variable('ib_torque_x', t_step, torque_x, nBodies) - call s_write_ib_variable('ib_torque_y', t_step, torque_y, nBodies) - call s_write_ib_variable('ib_torque_z', t_step, torque_z, nBodies) - call s_write_ib_variable('ib_vel_x', t_step, vel_x, nBodies) - call s_write_ib_variable('ib_vel_y', t_step, vel_y, nBodies) - call s_write_ib_variable('ib_vel_z', t_step, vel_z, nBodies) - call s_write_ib_variable('ib_omega_x', t_step, omega_x, nBodies) - call s_write_ib_variable('ib_omega_y', t_step, omega_y, nBodies) - call s_write_ib_variable('ib_omega_z', t_step, omega_z, nBodies) - call s_write_ib_variable('ib_angle_x', t_step, angle_x, nBodies) - call s_write_ib_variable('ib_angle_y', t_step, angle_y, nBodies) - call s_write_ib_variable('ib_angle_z', t_step, angle_z, nBodies) - call s_write_ib_variable('ib_diameter', t_step, ib_diameter, nBodies) - deallocate (ib_data, px, py, pz, force_x, force_y, force_z) deallocate (torque_x, torque_y, torque_z, vel_x, vel_y, vel_z) deallocate (omega_x, omega_y, omega_z, angle_x, angle_y, angle_z) @@ -1450,23 +1444,18 @@ contains !> Write a single IB point-variable to the Silo database slave and master files. subroutine s_write_ib_variable(varname, t_step, data, nBodies) - character(len=*), intent(in) :: varname - integer, intent(in) :: t_step - real(wp), dimension(:), intent(in) :: data - integer, intent(in) :: nBodies - character(len=4*name_len), dimension(num_procs) :: var_names - integer, dimension(num_procs) :: var_types - integer :: ierr, i - - if (proc_rank == 0) then - do i = 1, num_procs - write (var_names(i), '(A,I0,A,I0,A)') '../p', i - 1, '/', t_step, '.silo:' // trim(varname) - var_types(i) = DB_POINTVAR - end do - err = DBSET2DSTRLEN(len(var_names(1))) - err = DBPUTMVAR(dbroot, trim(varname), len_trim(varname), num_procs, var_names, len_trim(var_names), var_types, & - & DB_F77NULL, ierr) - end if + character(len=*), intent(in) :: varname + integer, intent(in) :: t_step + real(wp), dimension(:), intent(in) :: data + integer, intent(in) :: nBodies + character(len=4*name_len) :: var_name_entry + integer :: var_type_entry, ierr + + write (var_name_entry, '(A,I0,A)') '../p0/', t_step, '.silo:' // trim(varname) + var_type_entry = DB_POINTVAR + err = DBSET2DSTRLEN(len(var_name_entry)) + err = DBPUTMVAR(dbroot, trim(varname), len_trim(varname), 1, var_name_entry, len_trim(var_name_entry), var_type_entry, & + & DB_F77NULL, ierr) err = DBPUTPV1(dbfile, trim(varname), len_trim(varname), 'ib_bodies', 9, data, nBodies, DB_DOUBLE, DB_F77NULL, ierr) diff --git a/src/simulation/m_collisions.fpp b/src/simulation/m_collisions.fpp index fe860cb055..fdb9d14395 100644 --- a/src/simulation/m_collisions.fpp +++ b/src/simulation/m_collisions.fpp @@ -19,7 +19,8 @@ module m_collisions implicit none - private; public :: s_apply_collision_forces, s_initialize_collisions_module, s_finalize_collisions_module + private; public :: s_apply_collision_forces, s_initialize_collisions_module, s_finalize_collisions_module, & + & f_local_rank_owns_location, f_neighborhood_ranks_own_location ! overlap distances for computing collisions integer, allocatable, dimension(:,:) :: collision_lookup real(wp), allocatable, dimension(:,:) :: wall_overlap_distances @@ -61,7 +62,8 @@ contains ! get is distance used in the force calculation with each IB and each wall call s_detect_wall_collisions() - call s_detect_ib_collisions(ghost_points, ib_markers, num_gps, num_considered_collisions) + ! call s_detect_ib_collisions(ghost_points, ib_markers, num_gps, num_considered_collisions) + call s_detect_ib_collisions_n2(num_considered_collisions) select case (collision_model) case (1) ! soft sphere model @@ -97,6 +99,8 @@ contains encoded_pid2 = collision_lookup(i, 4) call s_decode_patch_periodicity(encoded_pid1, pid1, xp1, yp1, zp1) call s_decode_patch_periodicity(encoded_pid2, pid2, xp2, yp2, zp2) + ! call s_get_neighborhood_idx(pid1, pid1) ! global patch ID -> local index call s_get_neighborhood_idx(pid2, pid2) + if (pid1 <= 0 .or. pid2 <= 0) cycle centroid_1(1) = patch_ib(pid1)%x_centroid + real(xp1, wp)*(x_domain%end - x_domain%beg) centroid_1(2) = patch_ib(pid1)%y_centroid + real(yp1, wp)*(y_domain%end - y_domain%beg) @@ -113,7 +117,7 @@ contains overlap_distance = patch_ib(pid1)%radius + patch_ib(pid2)%radius - norm2(normal_vector) if (overlap_distance > 0._wp) then ! if the two patches are close enough to collide normal_vector = normal_vector/norm2(normal_vector) - if (f_local_rank_owns_collision(centroid_1)) then + if (f_local_rank_owns_location(centroid_1)) then ! compute constants of the collision effective_mass = 1.0_wp/((1.0_wp/patch_ib(pid1)%mass) + (1._wp/(patch_ib(pid2)%mass))) k = spring_stiffness*effective_mass @@ -192,7 +196,7 @@ contains ! ensure the local rank owns that collision before proceeding collision_location = [patch_ib(patch_id)%x_centroid, patch_ib(patch_id)%y_centroid, 0._wp] if (num_dims == 3) collision_location(3) = patch_ib(patch_id)%z_centroid - if (f_local_rank_owns_collision(collision_location)) then + if (f_local_rank_owns_location(collision_location)) then k = spring_stiffness*patch_ib(patch_id)%mass eta = damping_parameter*patch_ib(patch_id)%mass @@ -350,6 +354,8 @@ contains collision_lookup(current_collisions, 1) = pid1 collision_lookup(current_collisions, 2) = pid2 + collision_lookup(current_collisions, 3) = pid1 + collision_lookup(current_collisions, 4) = pid2 end if end do end do @@ -397,45 +403,78 @@ contains end subroutine s_detect_wall_collisions !> @brief function checks if this local MPI processor owns this specific collision - function f_local_rank_owns_collision(collision_location) result(owns_collision) + function f_local_rank_owns_location(location) result(owns_collision) $:GPU_ROUTINE(parallelism='[seq]') - real(wp), dimension(3), intent(in) :: collision_location + real(wp), dimension(3), intent(in) :: location logical :: owns_collision real(wp), dimension(3) :: projected_location - #:if defined('MFC_MPI') - if (num_procs == 1) then - owns_collision = .true. - else - projected_location(:) = collision_location(:) - - ! catch the edge case where th collision lies just outside the computational domain - #:for X, ID in [('x', 1), ('y', 2), ('z', 3)] - if (num_dims >= ${ID}$) then - if (ib_bc_${X}$%beg /= BC_PERIODIC) then - ! if it is outside the domain in one direction, project it somewhere inside so at least one rank owns it - if (collision_location(${ID}$) < ${X}$_domain%beg) then - projected_location(${ID}$) = ${X}$_domain%beg - else if (${X}$_domain%end < collision_location(${ID}$)) then - projected_location(${ID}$) = ${X}$_domain%end - 1.0e-10_wp - end if + owns_collision = .true. + +#ifdef MFC_MPI + if (num_procs > 1) then + projected_location(:) = location(:) + + ! catch the edge case where th collision lies just outside the computational domain + #:for X, ID, DIM in [('x', 1, 'm'), ('y', 2, 'n'), ('z', 3, 'p')] + if (num_dims >= ${ID}$) then + if (ib_bc_${X}$%beg /= BC_PERIODIC) then + ! if it is outside the domain in one direction, project it somewhere inside so at least one rank owns it + if (location(${ID}$) < ${X}$_domain%beg) then + projected_location(${ID}$) = ${X}$_domain%beg + else if (${X}$_domain%end < location(${ID}$)) then + projected_location(${ID}$) = ${X}$_domain%end - 1.0e-10_wp end if end if - #:endfor + owns_collision = owns_collision .and. ${X}$_cb(-1) <= projected_location(${ID}$) & + & .and. projected_location(${ID}$) < ${X}$_cb(${DIM}$) + end if + #:endfor + end if +#endif - ! the object that contains the collision location owns the collisions - owns_collision = x_cb(-1) <= projected_location(1) .and. projected_location(1) < x_cb(m) - owns_collision = owns_collision .and. y_cb(-1) <= projected_location(2) .and. projected_location(2) < y_cb(n) - if (num_dims == 3) owns_collision = owns_collision .and. z_cb(-1) <= projected_location(3) & - & .and. projected_location(3) < z_cb(p) - end if - #:else + end function f_local_rank_owns_location + + !> @brief function checks if this local MPI processor owns this specific collision + function f_neighborhood_ranks_own_location(location) result(owns_collision) + + $:GPU_ROUTINE(parallelism='[seq]') + + real(wp), dimension(3), intent(in) :: location + logical :: owns_collision, periodic_owner + real(wp) :: temp_neighbor_domain + integer :: i + + owns_collision = .true. + +#ifdef MFC_MPI + if (num_procs > 2) then + ! catch the edge case where th collision lies just outside the computational domain owns_collision = .true. - #:endif + #:for X, ID in [('x', 1), ('y', 2,), ('z', 3,)] + if (num_dims >= ${ID}$) then + if (ib_bc_${X}$%beg == BC_PERIODIC .and. neighbor_domain_${X}$%beg >= neighbor_domain_${X}$%end) then + ! project right side to the left + temp_neighbor_domain = neighbor_domain_${X}$%end + (${X}$_domain%end - ${X}$_domain%beg) + periodic_owner = neighbor_domain_${X}$%beg <= location(${ID}$) .and. location(${ID}$) < temp_neighbor_domain + ! project the left side to the right + temp_neighbor_domain = neighbor_domain_${X}$%beg - (${X}$_domain%end - ${X}$_domain%beg) + periodic_owner = periodic_owner .or. (temp_neighbor_domain <= location(${ID}$) .and. location(${ID}$) & + & < neighbor_domain_${X}$%end) + + owns_collision = owns_collision .and. periodic_owner + else + owns_collision = owns_collision .and. neighbor_domain_${X}$%beg <= location(${ID}$) .and. location(${ID}$) & + & < neighbor_domain_${X}$%end + end if + end if + #:endfor + end if +#endif - end function f_local_rank_owns_collision + end function f_neighborhood_ranks_own_location subroutine s_finalize_collisions_module() diff --git a/src/simulation/m_compute_levelset.fpp b/src/simulation/m_compute_levelset.fpp index 9e7519790c..6b4b31d7cb 100644 --- a/src/simulation/m_compute_levelset.fpp +++ b/src/simulation/m_compute_levelset.fpp @@ -88,8 +88,8 @@ contains radius = patch_ib(ib_patch_id)%radius - dist_vec(1) = x_cc(i) - patch_ib(ib_patch_id)%x_centroid - real(gp%x_periodicity, wp)*(x_domain%end - x_domain%beg) - dist_vec(2) = y_cc(j) - patch_ib(ib_patch_id)%y_centroid - real(gp%y_periodicity, wp)*(y_domain%end - y_domain%beg) + dist_vec(1) = x_cc(i) - (patch_ib(ib_patch_id)%x_centroid + real(gp%x_periodicity, wp)*(x_domain%end - x_domain%beg)) + dist_vec(2) = y_cc(j) - (patch_ib(ib_patch_id)%y_centroid + real(gp%y_periodicity, wp)*(y_domain%end - y_domain%beg)) dist_vec(3) = 0._wp dist = sqrt(sum(dist_vec**2)) diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index a5fdb3ef81..effc020325 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -909,7 +909,7 @@ contains integer :: ifile, ierr integer, dimension(MPI_STATUS_SIZE) :: status logical :: file_exist - integer :: i + integer :: i, ib_idx integer, parameter :: NFIELDS_PER_IB = 20 real(wp) :: ib_buf(NFIELDS_PER_IB) @@ -923,19 +923,6 @@ contains end if call s_mpi_barrier() - ! Divide num_ibs across num_procs - nibs_per_rank = num_ibs/num_procs - remainder = mod(num_ibs, num_procs) - - ! Ranks < remainder get one extra IB - if (proc_rank < remainder) then - ib_start = proc_rank*(nibs_per_rank + 1) + 1 - ib_end = ib_start + nibs_per_rank ! nibs_per_rank + 1 total - else - ib_start = remainder*(nibs_per_rank + 1) + (proc_rank - remainder)*nibs_per_rank + 1 - ib_end = ib_start + nibs_per_rank - 1 - end if - write (file_loc, '(A,I0,A)') '/restart_data/ib_state_', t_step, '.dat' file_loc = trim(case_dir) // trim(file_loc) @@ -947,20 +934,21 @@ contains call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, ior(MPI_MODE_WRONLY, MPI_MODE_CREATE), mpi_info_int, ifile, ierr) - do i = ib_start, ib_end + do i = 1, num_local_ibs + ib_idx = local_ib_patch_ids(i) ib_buf(1) = mytime - ib_buf(2:4) = patch_ib(i)%force(1:3) - ib_buf(5:7) = patch_ib(i)%torque(1:3) - ib_buf(8:10) = patch_ib(i)%vel(1:3) - ib_buf(11:13) = patch_ib(i)%angular_vel(1:3) - ib_buf(14:16) = patch_ib(i)%angles(1:3) - ib_buf(17) = patch_ib(i)%x_centroid - ib_buf(18) = patch_ib(i)%y_centroid - ib_buf(19) = patch_ib(i)%z_centroid - ib_buf(20) = patch_ib(i)%radius + ib_buf(2:4) = patch_ib(ib_idx)%force(1:3) + ib_buf(5:7) = patch_ib(ib_idx)%torque(1:3) + ib_buf(8:10) = patch_ib(ib_idx)%vel(1:3) + ib_buf(11:13) = patch_ib(ib_idx)%angular_vel(1:3) + ib_buf(14:16) = patch_ib(ib_idx)%angles(1:3) + ib_buf(17) = patch_ib(ib_idx)%x_centroid + ib_buf(18) = patch_ib(ib_idx)%y_centroid + ib_buf(19) = patch_ib(ib_idx)%z_centroid + ib_buf(20) = patch_ib(ib_idx)%radius ! Global IB index (i) determines position in file - disp = int(i - 1, MPI_OFFSET_KIND)*int(NFIELDS_PER_IB, MPI_OFFSET_KIND)*WP_MOK + disp = int(patch_ib(ib_idx)%gbl_patch_id - 1, MPI_OFFSET_KIND)*int(NFIELDS_PER_IB, MPI_OFFSET_KIND)*WP_MOK call MPI_FILE_WRITE_AT(ifile, disp, ib_buf, NFIELDS_PER_IB, mpi_p, status, ierr) end do diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index 8842759cbe..9ec0c2da5b 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -28,6 +28,7 @@ module m_global_parameters integer :: t_step_old !< Existing IC/grid folder ! Computational Domain Parameters integer :: proc_rank !< Rank of the local processor + $:GPU_DECLARE(create='[num_procs, proc_rank]') !> @name Number of cells in the x-, y- and z-directions, respectively !> @{ integer :: m, n, p @@ -232,7 +233,8 @@ module m_global_parameters $:GPU_DECLARE(create='[ib_bc_x, ib_bc_y, ib_bc_z]') #endif type(bounds_info) :: x_domain, y_domain, z_domain - $:GPU_DECLARE(create='[x_domain, y_domain, z_domain]') + type(bounds_info) :: neighbor_domain_x, neighbor_domain_y, neighbor_domain_z + $:GPU_DECLARE(create='[x_domain, y_domain, z_domain, neighbor_domain_x, neighbor_domain_y, neighbor_domain_z]') real(wp) :: x_a, y_a, z_a real(wp) :: x_b, y_b, z_b logical :: parallel_io !< Format of the data files @@ -339,18 +341,37 @@ module m_global_parameters !> @name Immersed Boundaries !> @{ - logical :: ib - integer :: num_ibs - integer :: collision_model - real(wp) :: coefficient_of_restitution - real(wp) :: collision_time - real(wp) :: ib_coefficient_of_friction - logical :: ib_state_wrt - type(ib_patch_parameters), dimension(num_ib_patches_max) :: patch_ib !< Immersed boundary patch parameters - type(vec3_dt), allocatable, dimension(:) :: airfoil_grid_u, airfoil_grid_l - integer :: Np - - $:GPU_DECLARE(create='[ib, num_ibs, patch_ib, Np, airfoil_grid_u, airfoil_grid_l]') + logical :: ib + integer :: num_ibs !< number of IBs that the current processor is aware of + integer :: num_gbl_ibs !< number of IBs in the overall simulation + integer :: num_local_ibs !< number of IBs that lie inside the processor domain + integer :: ib_neighborhood_radius !< neighborhood radius in ranks (1 = immediate neighbors) + integer :: collision_model + real(wp) :: coefficient_of_restitution + real(wp) :: collision_time + real(wp) :: ib_coefficient_of_friction + logical :: ib_state_wrt + type(ib_patch_parameters), allocatable, dimension(:) :: patch_ib !< Immersed boundary patch parameters + integer, dimension(num_local_ibs_max) :: local_ib_patch_ids !< lookup table of IBs in the local compute domain + + type particle_bed_parameters + real(wp) :: x_centroid, y_centroid, z_centroid !< Center of the particle bed region + real(wp) :: length_x, length_y, length_z !< Dimensions of the particle bed region + integer :: num_particles !< Number of particles to generate + real(wp) :: radius !< Particle radius + real(wp) :: mass !< Particle mass + real(wp) :: min_spacing !< Minimum surface-to-surface gap (particle centers are 2*radius + min_spacing apart) + integer :: moving_ibm !< Motion flag: 0=static, 1=moving (forces), 2=forced path + integer :: seed !< Random seed for reproducible placement + end type particle_bed_parameters + + integer :: num_particle_beds !< Number of particle bed specifications + type(particle_bed_parameters), dimension(num_particle_beds_max) :: particle_bed !< Particle bed specifications + integer, allocatable, dimension(:,:,:) :: ib_neighbor_ranks !< MPI ranks of neighborhood domains, indexed (-N:N,-N:N,-N:N) + type(vec3_dt), allocatable, dimension(:) :: airfoil_grid_u, airfoil_grid_l + integer :: Np + + $:GPU_DECLARE(create='[ib, num_ibs, num_gbl_ibs, num_local_ibs, patch_ib, Np, airfoil_grid_u, airfoil_grid_l, local_ib_patch_ids]') $:GPU_DECLARE(create='[ib_coefficient_of_friction]') !> @} @@ -638,6 +659,7 @@ contains ! Immersed Boundaries ib = .false. num_ibs = dflt_int + ib_neighborhood_radius = 1 collision_model = 0 coefficient_of_restitution = dflt_real collision_time = dflt_real @@ -791,7 +813,25 @@ contains relativity = .false. #:endif - do i = 1, num_ib_patches_max + num_particle_beds = 0 + do i = 1, num_particle_beds_max + particle_bed(i)%x_centroid = 0._wp + particle_bed(i)%y_centroid = 0._wp + particle_bed(i)%z_centroid = 0._wp + particle_bed(i)%length_x = dflt_real + particle_bed(i)%length_y = dflt_real + particle_bed(i)%length_z = dflt_real + particle_bed(i)%num_particles = 0 + particle_bed(i)%radius = dflt_real + particle_bed(i)%mass = dflt_real + particle_bed(i)%min_spacing = 0._wp + particle_bed(i)%moving_ibm = 0 + particle_bed(i)%seed = 0 + end do + + allocate (patch_ib(num_ib_patches_max_namelist)) + do i = 1, num_ib_patches_max_namelist + patch_ib(i)%gbl_patch_id = i patch_ib(i)%geometry = dflt_int patch_ib(i)%x_centroid = 0._wp patch_ib(i)%y_centroid = 0._wp diff --git a/src/simulation/m_ib_patches.fpp b/src/simulation/m_ib_patches.fpp index b21483e097..1b41bfc62c 100644 --- a/src/simulation/m_ib_patches.fpp +++ b/src/simulation/m_ib_patches.fpp @@ -104,7 +104,7 @@ contains radius = patch_ib(patch_id)%radius ! encode the periodicity information into the patch_id - call s_encode_patch_periodicity(patch_id, xp, yp, 0, encoded_patch_id) + call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, 0, encoded_patch_id) ! find the indices to the left and right of the IB in i, j, k il = -gp_layers - 1 @@ -216,7 +216,7 @@ contains end if ! encode the periodicity information into the patch_id - call s_encode_patch_periodicity(patch_id, xp, yp, 0, encoded_patch_id) + call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, 0, encoded_patch_id) ! find the indices to the left and right of the IB in i, j, k il = -gp_layers - 1 @@ -371,7 +371,7 @@ contains end if ! encode the periodicity information into the patch_id - call s_encode_patch_periodicity(patch_id, xp, yp, zp, encoded_patch_id) + call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, zp, encoded_patch_id) ! find the indices to the left and right of the IB in i, j, k il = -gp_layers - 1 @@ -461,7 +461,7 @@ contains inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:) ! encode the periodicity information into the patch_id - call s_encode_patch_periodicity(patch_id, xp, yp, 0, encoded_patch_id) + call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, 0, encoded_patch_id) ! find the indices to the left and right of the IB in i, j, k il = -gp_layers - 1 @@ -505,16 +505,12 @@ contains real(wp) :: radius real(wp), dimension(1:3) :: center - ! Variables to initialize the pressure field that corresponds to the bubble-collapse test case found in Tiwari et al. (2013) - - ! Transferring spherical patch's radius, centroid, smoothing patch identity and smoothing coefficient information - center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg) center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg) center(3) = patch_ib(patch_id)%z_centroid + real(zp, wp)*(z_domain%end - z_domain%beg) radius = patch_ib(patch_id)%radius - ! completely skip particles no in the domain + ! completely skip particles not in the domain if (center(1) - radius > x_cc(m + gp_layers + 1) .or. center(1) + radius < x_cc(-gp_layers - 1) .or. center(2) & & - radius > y_cc(n + gp_layers + 1) .or. center(2) + radius < y_cc(-gp_layers - 1) .or. center(3) - radius > z_cc(p & & + gp_layers + 1) .or. center(3) + radius < z_cc(-gp_layers - 1)) then @@ -522,7 +518,7 @@ contains end if ! encode the periodicity information into the patch_id - call s_encode_patch_periodicity(patch_id, xp, yp, zp, encoded_patch_id) + call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, zp, encoded_patch_id) ! find the indices to the left and right of the IB in i, j, k il = -gp_layers - 1 @@ -575,7 +571,7 @@ contains inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:) ! encode the periodicity information into the patch_id - call s_encode_patch_periodicity(patch_id, xp, yp, zp, encoded_patch_id) + call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, zp, encoded_patch_id) ! find the indices to the left and right of the IB in i, j, k il = -gp_layers - 1 @@ -638,7 +634,7 @@ contains inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:) ! encode the periodicity information into the patch_id - call s_encode_patch_periodicity(patch_id, xp, yp, zp, encoded_patch_id) + call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, zp, encoded_patch_id) il = -gp_layers - 1 jl = -gp_layers - 1 @@ -700,7 +696,7 @@ contains inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:) ! encode the periodicity information into the patch_id - call s_encode_patch_periodicity(patch_id, xp, yp, 0, encoded_patch_id) + call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, 0, encoded_patch_id) ! find the indices to the left and right of the IB in i, j, k il = -gp_layers - 1 @@ -738,7 +734,7 @@ contains integer, intent(in) :: xp, yp !< integers containing the periodicity projection information integer :: i, j, il, ir, jl, jr !< Generic loop iterators integer :: spc, encoded_patch_id - integer :: cx, cy + integer :: cx, cy, gbl_patch_id real(wp) :: lx(2), ly(2) real(wp), dimension(1:2) :: bbox_min, bbox_max real(wp), dimension(1:3) :: local_corner, world_corner @@ -747,6 +743,7 @@ contains real(wp), dimension(1:3) :: center, xy_local real(wp), dimension(1:3,1:3) :: inverse_rotation, rotation + gbl_patch_id = patch_ib(patch_id)%gbl_patch_id center = 0._wp center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg) center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg) @@ -757,7 +754,7 @@ contains threshold = patch_ib(patch_id)%model_threshold ! encode the periodicity information into the patch_id - call s_encode_patch_periodicity(patch_id, xp, yp, 0, encoded_patch_id) + call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, 0, encoded_patch_id) il = -gp_layers - 1 jl = -gp_layers - 1 @@ -765,10 +762,10 @@ contains jr = n + gp_layers + 1 ! Local-space bounding box extents (min=1, max=2 in the third index) - lx(1) = stl_bounding_boxes(patch_id, 1, 1) + offset(1) - lx(2) = stl_bounding_boxes(patch_id, 1, 3) + offset(1) - ly(1) = stl_bounding_boxes(patch_id, 2, 1) + offset(2) - ly(2) = stl_bounding_boxes(patch_id, 2, 3) + offset(2) + lx(1) = stl_bounding_boxes(gbl_patch_id, 1, 1) + offset(1) + lx(2) = stl_bounding_boxes(gbl_patch_id, 1, 3) + offset(1) + ly(1) = stl_bounding_boxes(gbl_patch_id, 2, 1) + offset(2) + ly(2) = stl_bounding_boxes(gbl_patch_id, 2, 3) + offset(2) bbox_min = 1e12 bbox_max = -1e12 @@ -795,7 +792,7 @@ contains xy_local = matmul(inverse_rotation, xy_local) xy_local = xy_local - offset - eta = f_model_is_inside_flat(gpu_ntrs(patch_id), patch_id, xy_local) + eta = f_model_is_inside_flat(gpu_ntrs(gbl_patch_id), gbl_patch_id, xy_local) ! Reading STL boundary vertices and compute the levelset and levelset_norm if (eta > threshold) then @@ -814,7 +811,7 @@ contains type(integer_field), intent(inout) :: ib_markers integer, intent(in) :: xp, yp, zp !< integers containing the periodicity projection information integer :: i, j, k, il, ir, jl, jr, kl, kr !< Generic loop iterators - integer :: spc, encoded_patch_id + integer :: spc, encoded_patch_id, gbl_patch_id real(wp) :: eta, threshold real(wp), dimension(1:3) :: offset real(wp), dimension(1:3) :: center, xyz_local @@ -823,6 +820,7 @@ contains real(wp) :: lx(2), ly(2), lz(2) real(wp), dimension(1:3) :: bbox_min, bbox_max, local_corner, world_corner + gbl_patch_id = patch_ib(patch_id)%gbl_patch_id center = 0._wp center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg) center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg) @@ -834,7 +832,7 @@ contains rotation(:,:) = patch_ib(patch_id)%rotation_matrix(:,:) ! encode the periodicity information into the patch_id - call s_encode_patch_periodicity(patch_id, xp, yp, zp, encoded_patch_id) + call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, zp, encoded_patch_id) il = -gp_layers - 1 jl = -gp_layers - 1 @@ -844,12 +842,12 @@ contains kr = p + gp_layers + 1 ! Local-space bounding box extents (min=1, max=2 in the third index) - lx(1) = stl_bounding_boxes(patch_id, 1, 1) + offset(1) - lx(2) = stl_bounding_boxes(patch_id, 1, 3) + offset(1) - ly(1) = stl_bounding_boxes(patch_id, 2, 1) + offset(2) - ly(2) = stl_bounding_boxes(patch_id, 2, 3) + offset(2) - lz(1) = stl_bounding_boxes(patch_id, 3, 1) + offset(3) - lz(2) = stl_bounding_boxes(patch_id, 3, 3) + offset(3) + lx(1) = stl_bounding_boxes(gbl_patch_id, 1, 1) + offset(1) + lx(2) = stl_bounding_boxes(gbl_patch_id, 1, 3) + offset(1) + ly(1) = stl_bounding_boxes(gbl_patch_id, 2, 1) + offset(2) + ly(2) = stl_bounding_boxes(gbl_patch_id, 2, 3) + offset(2) + lz(1) = stl_bounding_boxes(gbl_patch_id, 3, 1) + offset(3) + lz(2) = stl_bounding_boxes(gbl_patch_id, 3, 3) + offset(3) bbox_min = 1e12 bbox_max = -1e12 @@ -882,7 +880,7 @@ contains xyz_local = matmul(inverse_rotation, xyz_local) xyz_local = xyz_local - offset - eta = f_model_is_inside_flat(gpu_ntrs(patch_id), patch_id, xyz_local) + eta = f_model_is_inside_flat(gpu_ntrs(gbl_patch_id), gbl_patch_id, xyz_local) if (eta > patch_ib(patch_id)%model_threshold) then ib_markers%sf(i, j, k) = encoded_patch_id @@ -991,7 +989,7 @@ contains temp_y_per = y_periodicity; if (y_periodicity == -1) temp_y_per = 2 temp_z_per = z_periodicity; if (z_periodicity == -1) temp_z_per = 2 - offset = (num_ibs + 1)*temp_x_per + 3*(num_ibs + 1)*temp_y_per + 9*(num_ibs + 1)*temp_z_per + offset = (num_gbl_ibs + 1)*temp_x_per + 3*(num_gbl_ibs + 1)*temp_y_per + 9*(num_gbl_ibs + 1)*temp_z_per encoded_patch_id = patch_id + offset end subroutine s_encode_patch_periodicity @@ -1006,7 +1004,7 @@ contains integer, intent(out), optional :: x_periodicity, y_periodicity, z_periodicity integer :: offset, remainder, xp, yp, zp, base - base = num_ibs + 1 + base = num_gbl_ibs + 1 patch_id = mod(encoded_patch_id - 1, base) + 1 offset = (encoded_patch_id - patch_id)/base diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 39ec6b1470..6457e943dc 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -25,7 +25,7 @@ module m_ibm private :: s_compute_image_points, s_compute_interpolation_coeffs, s_interpolate_image_point, s_find_ghost_points, & & s_find_num_ghost_points - ; public :: s_initialize_ibm_module, s_ibm_setup, s_ibm_correct_state, s_finalize_ibm_module + ; public :: ib_gbl_idx_lookup, s_initialize_ibm_module, s_ibm_setup, s_ibm_correct_state, s_finalize_ibm_module type(integer_field), public :: ib_markers $:GPU_DECLARE(create='[ib_markers]') @@ -33,6 +33,9 @@ module m_ibm type(ghost_point), dimension(:), allocatable :: ghost_points $:GPU_DECLARE(create='[ghost_points]') + integer, dimension(:), allocatable :: ib_gbl_idx_lookup + $:GPU_DECLARE(create='[ib_gbl_idx_lookup]') + integer :: num_gps !< Number of ghost points #if defined(MFC_OpenACC) $:GPU_DECLARE(create='[gp_layers, num_gps]') @@ -41,6 +44,11 @@ module m_ibm #endif logical :: moving_immersed_boundary_flag + ! IB MPI buffers + integer, allocatable :: send_ids(:), recv_ids(:) + real(wp), allocatable :: send_ft(:,:), recv_ft(:,:) + real(wp), allocatable :: recv_forces_snap(:,:), recv_torques_snap(:,:) + contains !> Allocates memory for the variables in the IBM module @@ -71,26 +79,31 @@ contains call nvtxStartRange("SETUP-IBM-MODULE") ! do all set up for moving immersed boundaries - moving_immersed_boundary_flag = .false. do i = 1, num_ibs if (patch_ib(i)%moving_ibm /= 0) then call s_compute_moment_of_inertia(i, patch_ib(i)%angular_vel) - moving_immersed_boundary_flag = .true. end if call s_update_ib_rotation_matrix(i) end do $:GPU_UPDATE(device='[patch_ib(1:num_ibs)]') + ! allocate some arrays for MPI communication, if required by this simulation +#ifdef MFC_MPI + if (num_procs > 1) then + @:ALLOCATE(send_ids(size(patch_ib)), send_ft(6, size(patch_ib))) + allocate (recv_forces_snap(size(patch_ib), 3), recv_torques_snap(size(patch_ib), 3), recv_ids(size(patch_ib)), & + & recv_ft(6, size(patch_ib))) + end if +#endif + ! GPU routines require updated cell centers - $:GPU_UPDATE(device='[num_ibs, x_cc, y_cc, dx, dy, x_domain, y_domain, ib_bc_x%beg, ib_bc_y%beg]') + $:GPU_UPDATE(device='[num_ibs, num_gbl_ibs, x_cc, y_cc, dx, dy, x_domain, y_domain, ib_bc_x%beg, ib_bc_y%beg]') if (p /= 0) then $:GPU_UPDATE(device='[z_cc, dz, z_domain, ib_bc_z%beg]') end if + call s_update_ib_lookup() - ! allocate STL models - call s_instantiate_STL_models() - - ! recompute the new ib_patch locations and broadcast them. + ! recompute the new ib_patch locations ib_markers%sf = 0._wp $:GPU_UPDATE(device='[ib_markers%sf]') call s_apply_ib_patches(ib_markers) @@ -174,17 +187,21 @@ contains do j = 0, m patch_id = ib_markers%sf(j, k, l) if (patch_id /= 0) then - q_prim_vf(eqn_idx%E)%sf(j, k, l) = 1._wp - rho = 0._wp - do i = 1, num_fluids - rho = rho + q_prim_vf(eqn_idx%cont%beg + i - 1)%sf(j, k, l) - end do + call s_decode_patch_periodicity(patch_id, patch_id) + call s_get_neighborhood_idx(patch_id, patch_id) + if (patch_id > 0) then + q_prim_vf(eqn_idx%E)%sf(j, k, l) = 1._wp + rho = 0._wp + do i = 1, num_fluids + rho = rho + q_prim_vf(eqn_idx%cont%beg + i - 1)%sf(j, k, l) + end do - ! Sets the momentum - do i = 1, num_dims - q_cons_vf(eqn_idx%mom%beg + i - 1)%sf(j, k, l) = patch_ib(patch_id)%vel(i)*rho - q_prim_vf(eqn_idx%mom%beg + i - 1)%sf(j, k, l) = patch_ib(patch_id)%vel(i) - end do + ! Sets the momentum + do i = 1, num_dims + q_cons_vf(eqn_idx%mom%beg + i - 1)%sf(j, k, l) = patch_ib(patch_id)%vel(i)*rho + q_prim_vf(eqn_idx%mom%beg + i - 1)%sf(j, k, l) = patch_ib(patch_id)%vel(i) + end do + end if ! patch_id > 0 end if end do end do @@ -479,7 +496,7 @@ contains end do $:END_GPU_PARALLEL_LOOP() - if (bounds_error) error stop "Ghost Point and Image Point on Different Processors. Exiting" + @:PROHIBIT(bounds_error, "Ghost Point and Image Point on Different Processors. Exiting") end subroutine s_compute_image_points @@ -535,7 +552,7 @@ contains integer :: i, j, k, ii, jj, kk, gp_layers_z !< Iterator variables integer :: xp, yp, zp !< periodicities integer :: count, count_i, local_idx - integer :: patch_id, encoded_patch_id + integer :: patch_id, encoded_patch_id, neighborhood_patch_id logical :: is_gp count = 0 @@ -543,8 +560,9 @@ contains gp_layers_z = gp_layers if (p == 0) gp_layers_z = 0 - $:GPU_PARALLEL_LOOP(private='[i, j, k, ii, jj, kk, is_gp, local_idx, patch_id, encoded_patch_id, xp, yp, zp]', & - & copyin='[count, count_i, x_domain, y_domain, z_domain]', firstprivate='[gp_layers, gp_layers_z]', collapse=3) + $:GPU_PARALLEL_LOOP(private='[i, j, k, ii, jj, kk, is_gp, local_idx, patch_id, encoded_patch_id, neighborhood_patch_id, & + & xp, yp, zp]', copyin='[count, count_i, x_domain, y_domain, z_domain]', firstprivate='[gp_layers, & + & gp_layers_z]', collapse=3) do i = 0, m do j = 0, n do k = 0, p @@ -571,11 +589,12 @@ contains ghost_points_in(local_idx)%loc = [i, j, k] encoded_patch_id = ib_markers%sf(i, j, k) call s_decode_patch_periodicity(encoded_patch_id, patch_id, xp, yp, zp) - ghost_points_in(local_idx)%ib_patch_id = patch_id + call s_get_neighborhood_idx(patch_id, neighborhood_patch_id) + ghost_points_in(local_idx)%ib_patch_id = neighborhood_patch_id ghost_points_in(local_idx)%x_periodicity = xp ghost_points_in(local_idx)%y_periodicity = yp ghost_points_in(local_idx)%z_periodicity = zp - ghost_points_in(local_idx)%slip = patch_ib(patch_id)%slip + ghost_points_in(local_idx)%slip = patch_ib(neighborhood_patch_id)%slip if ((x_cc(i) - dx(i)) < x_domain%beg) then ghost_points_in(local_idx)%DB(1) = -1 @@ -718,9 +737,10 @@ contains !> Interpolate primitive variables to a ghost point's image point using bilinear or trilinear interpolation subroutine s_interpolate_image_point(q_prim_vf, gp, alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP, r_IP, v_IP, pb_IP, mv_IP, & + & nmom_IP, pb_in, mv_in, presb_IP, massv_IP) - & nmom_IP, pb_in, mv_in, presb_IP, massv_IP) $:GPU_ROUTINE(parallelism='[seq]') + type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf !< Primitive Variables real(stp), optional, dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), intent(in) :: pb_in, mv_in type(ghost_point), intent(in) :: gp @@ -923,83 +943,91 @@ contains encoded_ib_idx = ib_markers%sf(i, j, k) if (encoded_ib_idx /= 0) then call s_decode_patch_periodicity(encoded_ib_idx, ib_idx) - - ! get the vector pointing to the grid cell from the IB centroid - if (num_dims == 3) then - radial_vector = [x_cc(i), y_cc(j), z_cc(k)] - [patch_ib(ib_idx)%x_centroid, & - & patch_ib(ib_idx)%y_centroid, patch_ib(ib_idx)%z_centroid] - else - radial_vector = [x_cc(i), y_cc(j), 0._wp] - [patch_ib(ib_idx)%x_centroid, & - & patch_ib(ib_idx)%y_centroid, 0._wp] - end if - dx = x_cc(i + 1) - x_cc(i) - dy = y_cc(j + 1) - y_cc(j) - - local_force_contribution(:) = 0._wp - do fluid_idx = 0, num_fluids - 1 - ! Get the pressure contribution to force via a finite difference to compute the 2D components of the - ! gradient of the pressure and cell volume - local_force_contribution(1) = local_force_contribution(1) - (q_prim_vf(eqn_idx%E + fluid_idx)%sf(i & - & + 1, j, k) - q_prim_vf(eqn_idx%E + fluid_idx)%sf(i - 1, j, k))/(2._wp*dx) ! force is the negative pressure gradient - local_force_contribution(2) = local_force_contribution(2) - (q_prim_vf(eqn_idx%E + fluid_idx)%sf(i, & - & j + 1, k) - q_prim_vf(eqn_idx%E + fluid_idx)%sf(i, j - 1, k))/(2._wp*dy) - cell_volume = abs(dx*dy) - ! add the 3D component of the pressure gradient, if we are working in 3 dimensions + call s_get_neighborhood_idx(ib_idx, ib_idx) ! global patch ID -> local index + if (ib_idx > 0) then + ! get the vector pointing to the grid cell from the IB centroid if (num_dims == 3) then - dz = z_cc(k + 1) - z_cc(k) - local_force_contribution(3) = local_force_contribution(3) - (q_prim_vf(eqn_idx%E & - & + fluid_idx)%sf(i, j, k + 1) - q_prim_vf(eqn_idx%E + fluid_idx)%sf(i, & - & j, k - 1))/(2._wp*dz) - cell_volume = abs(cell_volume*dz) + radial_vector = [x_cc(i), y_cc(j), z_cc(k)] - [patch_ib(ib_idx)%x_centroid, & + & patch_ib(ib_idx)%y_centroid, patch_ib(ib_idx)%z_centroid] + else + radial_vector = [x_cc(i), y_cc(j), 0._wp] - [patch_ib(ib_idx)%x_centroid, & + & patch_ib(ib_idx)%y_centroid, 0._wp] end if - end do - - ! get the viscous stress and add its contribution if that is considered - if (viscous) then - ! compute the volume-weighted local dynamic viscosity - dynamic_viscosity = 0._wp - do fluid_idx = 1, num_fluids - ! local dynamic viscosity is the dynamic viscosity of the fluid times alpha of the fluid - dynamic_viscosity = dynamic_viscosity + (q_prim_vf(fluid_idx + eqn_idx%adv%beg - 1)%sf(i, j, & - & k)*dynamic_viscosities(fluid_idx)) + dx = x_cc(i + 1) - x_cc(i) + dy = y_cc(j + 1) - y_cc(j) + + local_force_contribution(:) = 0._wp + do fluid_idx = 0, num_fluids - 1 + ! Get the pressure contribution to force via a finite difference to compute the 2D components of the + ! gradient of the pressure and cell volume + local_force_contribution(1) = local_force_contribution(1) - (q_prim_vf(eqn_idx%E & + & + fluid_idx)%sf(i + 1, j, & + & k) - q_prim_vf(eqn_idx%E + fluid_idx)%sf(i - 1, j, k))/(2._wp*dx) ! force is the negative pressure gradient + local_force_contribution(2) = local_force_contribution(2) - (q_prim_vf(eqn_idx%E & + & + fluid_idx)%sf(i, j + 1, k) - q_prim_vf(eqn_idx%E + fluid_idx)%sf(i, & + & j - 1, k))/(2._wp*dy) + cell_volume = abs(dx*dy) + ! add the 3D component of the pressure gradient, if we are working in 3 dimensions + if (num_dims == 3) then + dz = z_cc(k + 1) - z_cc(k) + local_force_contribution(3) = local_force_contribution(3) - (q_prim_vf(eqn_idx%E & + & + fluid_idx)%sf(i, j, & + & k + 1) - q_prim_vf(eqn_idx%E + fluid_idx)%sf(i, j, k - 1))/(2._wp*dz) + cell_volume = abs(cell_volume*dz) + end if end do - ! get the linear force components first - call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i - 1, j, k) - call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i + 1, j, k) - ! get x derivative of the first-row of viscous stress tensor - viscous_stress_div(1,1:3) = (viscous_stress_div_2(1,1:3) - viscous_stress_div_1(1,1:3))/(2._wp*dx) - ! add the x components of the divergence to the force - local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(1,1:3) - - call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i, j - 1, k) - call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i, j + 1, k) - ! get y derivative of the second-row of viscous stress tensor - viscous_stress_div(2,1:3) = (viscous_stress_div_2(2,1:3) - viscous_stress_div_1(2,1:3))/(2._wp*dy) - ! add the y components of the divergence to the force - local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(2,1:3) + ! get the viscous stress and add its contribution if that is considered + if (viscous) then + ! compute the volume-weighted local dynamic viscosity + dynamic_viscosity = 0._wp + do fluid_idx = 1, num_fluids + ! local dynamic viscosity is the dynamic viscosity of the fluid times alpha of the fluid + dynamic_viscosity = dynamic_viscosity + (q_prim_vf(fluid_idx + eqn_idx%adv%beg - 1)%sf(i, j, & + & k)*dynamic_viscosities(fluid_idx)) + end do - if (num_dims == 3) then - call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i, j, & - & k - 1) - call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i, j, & - & k + 1) - ! get z derivative of the third-row of viscous stress tensor - viscous_stress_div(3,1:3) = (viscous_stress_div_2(3,1:3) - viscous_stress_div_1(3,1:3))/(2._wp*dz) - ! add the z components of the divergence to the force - local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(3,1:3) + ! get the linear force components first + call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i - 1, & + & j, k) + call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i + 1, & + & j, k) + ! get x derivative of the first-row of viscous stress tensor + viscous_stress_div(1,1:3) = (viscous_stress_div_2(1,1:3) - viscous_stress_div_1(1,1:3))/(2._wp*dx) + ! add the x components of the divergence to the force + local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(1,1:3) + + call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i, & + & j - 1, k) + call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i, & + & j + 1, k) + ! get y derivative of the second-row of viscous stress tensor + viscous_stress_div(2,1:3) = (viscous_stress_div_2(2,1:3) - viscous_stress_div_1(2,1:3))/(2._wp*dy) + ! add the y components of the divergence to the force + local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(2,1:3) + + if (num_dims == 3) then + call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i, & + & j, k - 1) + call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i, & + & j, k + 1) + viscous_stress_div(3,1:3) = (viscous_stress_div_2(3,1:3) - viscous_stress_div_1(3, & + & 1:3))/(2._wp*dz) + ! add the z components of the divergence to the force + local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(3,1:3) + end if end if - end if - call s_cross_product(radial_vector, local_force_contribution, local_torque_contribution) + call s_cross_product(radial_vector, local_force_contribution, local_torque_contribution) - ! Update the force and torque values atomically to prevent race conditions - do l = 1, 3 - $:GPU_ATOMIC(atomic='update') - forces(ib_idx, l) = forces(ib_idx, l) + (local_force_contribution(l)*cell_volume) - $:GPU_ATOMIC(atomic='update') - torques(ib_idx, l) = torques(ib_idx, l) + local_torque_contribution(l)*cell_volume - end do + ! Update the force and torque values atomically to prevent race conditions + do l = 1, 3 + $:GPU_ATOMIC(atomic='update') + forces(ib_idx, l) = forces(ib_idx, l) + (local_force_contribution(l)*cell_volume) + $:GPU_ATOMIC(atomic='update') + torques(ib_idx, l) = torques(ib_idx, l) + local_torque_contribution(l)*cell_volume + end do + end if ! ib_idx > 0 end if end do end do @@ -1008,9 +1036,8 @@ contains call s_apply_collision_forces(ghost_points, num_gps, ib_markers, forces, torques) - ! reduce the forces across all MPI ranks - call s_mpi_allreduce_vectors_sum(forces, forces, num_ibs, 3) - call s_mpi_allreduce_vectors_sum(torques, torques, num_ibs, 3) + ! reduce the forces across local neighborhood ranks + call s_communicate_ib_forces(forces, torques) ! consider body forces after reducing to avoid double counting do i = 1, num_ibs @@ -1035,25 +1062,12 @@ contains end subroutine s_compute_ib_forces - !> Finalize the IBM module - impure subroutine s_finalize_ibm_module() - - @:DEALLOCATE(ib_markers%sf) - if (allocated(airfoil_grid_u)) then - @:DEALLOCATE(airfoil_grid_u) - @:DEALLOCATE(airfoil_grid_l) - end if - - if (collision_model > 0) call s_finalize_collisions_module() - - end subroutine s_finalize_ibm_module - !> Computes the center of mass for IB patch types where we are unable to determine their center of mass analytically. !> These patches include things like NACA airfoils and STL models subroutine s_compute_centroid_offset(ib_marker) integer, intent(in) :: ib_marker - integer :: i, j, k, num_cells, num_cells_local + integer :: i, j, k, num_cells, num_cells_local, decoded_gbl_id real(wp), dimension(1:3) :: center_of_mass, center_of_mass_local ! Offset only needs to be computes for specific geometries @@ -1067,10 +1081,13 @@ contains do i = 0, m do j = 0, n do k = 0, p - if (ib_markers%sf(i, j, k) == ib_marker) then - num_cells_local = num_cells_local + 1 - center_of_mass_local = center_of_mass_local + [x_cc(i), y_cc(j), 0._wp] - if (num_dims == 3) center_of_mass_local(3) = center_of_mass_local(3) + z_cc(k) + if (ib_markers%sf(i, j, k) /= 0) then + call s_decode_patch_periodicity(ib_markers%sf(i, j, k), decoded_gbl_id) + if (decoded_gbl_id == patch_ib(ib_marker)%gbl_patch_id) then + num_cells_local = num_cells_local + 1 + center_of_mass_local = center_of_mass_local + [x_cc(i), y_cc(j), 0._wp] + if (num_dims == 3) center_of_mass_local(3) = center_of_mass_local(3) + z_cc(k) + end if end if end do end do @@ -1106,35 +1123,33 @@ contains end subroutine s_compute_centroid_offset !> Computes the moment of inertia for an immersed boundary - subroutine s_compute_moment_of_inertia(ib_marker, axis) + subroutine s_compute_moment_of_inertia(ib_idx, axis) real(wp), dimension(3), intent(in) :: axis !< the axis about which we compute the moment. Only required in 3D. - integer, intent(in) :: ib_marker + integer, intent(in) :: ib_idx real(wp) :: moment, distance_to_axis, cell_volume real(wp), dimension(3) :: position, closest_point_along_axis, vector_to_axis, normal_axis - integer :: i, j, k, count + integer :: i, j, k, count, ib_marker if (p == 0) then normal_axis = [0, 0, 1] else if (sqrt(sum(axis**2)) < sgm_eps) then ! if the object is not actually rotating at this time, return a dummy value and exit - patch_ib(ib_marker)%moment = 1._wp + patch_ib(ib_idx)%moment = 1._wp return else normal_axis = axis/sqrt(sum(axis)) end if ! if the IB is in 2D or a 3D sphere, we can compute this exactly - if (patch_ib(ib_marker)%geometry == 2) then ! circle - patch_ib(ib_marker)%moment = 0.5_wp*patch_ib(ib_marker)%mass*(patch_ib(ib_marker)%radius)**2 - else if (patch_ib(ib_marker)%geometry == 3) then ! rectangle - patch_ib(ib_marker)%moment = patch_ib(ib_marker)%mass*(patch_ib(ib_marker)%length_x**2 & - & + patch_ib(ib_marker)%length_y**2)/6._wp - else if (patch_ib(ib_marker)%geometry == 6) then ! ellipse - patch_ib(ib_marker)%moment = 0.0625_wp*patch_ib(ib_marker)%mass*(patch_ib(ib_marker)%length_x**2 & - & + patch_ib(ib_marker)%length_y**2) - else if (patch_ib(ib_marker)%geometry == 8) then ! sphere - patch_ib(ib_marker)%moment = 0.4*patch_ib(ib_marker)%mass*(patch_ib(ib_marker)%radius)**2 + if (patch_ib(ib_idx)%geometry == 2) then ! circle + patch_ib(ib_idx)%moment = 0.5_wp*patch_ib(ib_idx)%mass*(patch_ib(ib_idx)%radius)**2 + else if (patch_ib(ib_idx)%geometry == 3) then ! rectangle + patch_ib(ib_idx)%moment = patch_ib(ib_idx)%mass*(patch_ib(ib_idx)%length_x**2 + patch_ib(ib_idx) %length_y**2)/6._wp + else if (patch_ib(ib_idx)%geometry == 6) then ! ellipse + patch_ib(ib_idx)%moment = 0.0625_wp*patch_ib(ib_idx)%mass*(patch_ib(ib_idx)%length_x**2 + patch_ib(ib_idx) %length_y**2) + else if (patch_ib(ib_idx)%geometry == 8) then ! sphere + patch_ib(ib_idx)%moment = 0.4*patch_ib(ib_idx)%mass*(patch_ib(ib_idx)%radius)**2 else ! we do not have an analytic moment of inertia calculation and need to approximate it directly via a sum count = 0 moment = 0._wp @@ -1144,6 +1159,8 @@ contains cell_volume = cell_volume*(z_cc(1) - z_cc(0)) end if + ib_marker = patch_ib(ib_idx)%gbl_patch_id + $:GPU_PARALLEL_LOOP(private='[position, closest_point_along_axis, vector_to_axis, distance_to_axis]', copy='[moment, & & count]', copyin='[ib_marker, cell_volume, normal_axis]', collapse=3) do i = 0, m @@ -1154,12 +1171,12 @@ contains count = count + 1 ! increment the count of total cells in the boundary ! get the position in local coordinates so that the axis passes through 0, 0, 0 - if (p == 0) then - position = [x_cc(i), y_cc(j), 0._wp] - [patch_ib(ib_marker)%x_centroid, & - & patch_ib(ib_marker)%y_centroid, 0._wp] + if (num_dims < 3) then + position = [x_cc(i), y_cc(j), 0._wp] - [patch_ib(ib_idx)%x_centroid, patch_ib(ib_idx)%y_centroid, & + & 0._wp] else - position = [x_cc(i), y_cc(j), z_cc(k)] - [patch_ib(ib_marker)%x_centroid, & - & patch_ib(ib_marker)%y_centroid, patch_ib(ib_marker)%z_centroid] + position = [x_cc(i), y_cc(j), z_cc(k)] - [patch_ib(ib_idx)%x_centroid, & + & patch_ib(ib_idx)%y_centroid, patch_ib(ib_idx)%z_centroid] end if ! project the position along the axis to find the closest distance to the rotation axis @@ -1177,8 +1194,7 @@ contains $:END_GPU_PARALLEL_LOOP() ! write the final moment assuming the points are all uniform density - patch_ib(ib_marker)%moment = moment*patch_ib(ib_marker)%mass/(count*cell_volume) - $:GPU_UPDATE(device='[patch_ib(ib_marker)%moment]') + patch_ib(ib_idx)%moment = moment*patch_ib(ib_idx)%mass/(count*cell_volume) end if end subroutine s_compute_moment_of_inertia @@ -1223,4 +1239,319 @@ contains end subroutine s_wrap_periodic_ibs + !> @brief reasseses ownership of IBs and passes ownership of IBs to neighbor processors + !> Reduces forces and torques across the local neighborhood without a global allreduce. Accumulation phase: 2 passes per + !! dimension receiving from the low-index (-X) neighbor. Pass 1: add received values; save what was received as recv_snap. Pass + !! 2: send current (post-pass-1) values; add received; subtract recv_snap to remove double-counting of the direct contribution + !! already added in pass 1. Back-propagation phase: 2 passes per dimension receiving from the high-index (+X) neighbor, each + !! overwriting local forces with the neighbor's accumulated total. + subroutine s_communicate_ib_forces(forces, torques) + + real(wp), dimension(num_ibs, 3), intent(inout) :: forces, torques + +#ifdef MFC_MPI + integer :: i, j, k, pack_pos, unpack_pos, buf_size, ierr + integer :: send_neighbor, recv_neighbor, recv_count, tag + character(len=1), allocatable :: ib_force_send_buf(:), ib_force_recv_buf(:) + + if (num_procs == 1) return + + buf_size = storage_size(0)/8 + (storage_size(0)/8 + 6*storage_size(0._wp)/8)*size(patch_ib) + allocate (ib_force_send_buf(buf_size), ib_force_recv_buf(buf_size)) + + ! Accumulation phase: propagate contributions toward the high-index corner. + #:for X, ID in [('x', 1), ('y', 2), ('z', 3)] + if (num_dims >= ${ID}$) then + send_neighbor = merge(bc_${X}$%end, MPI_PROC_NULL, bc_${X}$%end >= 0) + recv_neighbor = merge(bc_${X}$%beg, MPI_PROC_NULL, bc_${X}$%beg >= 0) + + recv_forces_snap = 0._wp + recv_torques_snap = 0._wp + tag = 300 + + do k = 1, 2*ib_neighborhood_radius + ! send forces to +${X}$ neighbor; receive from -${X}$ neighbor. Add received values then + pack_pos = 0 + $:GPU_PARALLEL_LOOP(private='[i]', copyin='[forces, torques]') + do i = 1, num_ibs + send_ids(i) = patch_ib(i)%gbl_patch_id + send_ft(1:3,i) = forces(i,:) + send_ft(4:6,i) = torques(i,:) + end do + $:END_GPU_PARALLEL_LOOP() + $:GPU_UPDATE(host='[send_ids, send_ft]') + call MPI_PACK(num_ibs, 1, MPI_INTEGER, ib_force_send_buf, buf_size, pack_pos, MPI_COMM_WORLD, ierr) + call MPI_PACK(send_ids, num_ibs, MPI_INTEGER, ib_force_send_buf, buf_size, pack_pos, MPI_COMM_WORLD, ierr) + call MPI_PACK(send_ft, 6*num_ibs, mpi_p, ib_force_send_buf, buf_size, pack_pos, MPI_COMM_WORLD, ierr) + call MPI_SENDRECV(ib_force_send_buf, pack_pos, MPI_PACKED, send_neighbor, tag, ib_force_recv_buf, buf_size, & + & MPI_PACKED, recv_neighbor, tag, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) + + if (recv_neighbor /= MPI_PROC_NULL) then + unpack_pos = 0 + call MPI_UNPACK(ib_force_recv_buf, buf_size, unpack_pos, recv_count, 1, MPI_INTEGER, MPI_COMM_WORLD, ierr) + call MPI_UNPACK(ib_force_recv_buf, buf_size, unpack_pos, recv_ids, recv_count, MPI_INTEGER, & + & MPI_COMM_WORLD, ierr) + call MPI_UNPACK(ib_force_recv_buf, buf_size, unpack_pos, recv_ft, 6*recv_count, mpi_p, MPI_COMM_WORLD, ierr) + $:GPU_PARALLEL_LOOP(private='[i, j]', copyin='[recv_ft, recv_ids]', copy='[forces, torques, & + & recv_forces_snap, recv_torques_snap]') + do i = 1, recv_count + call s_get_neighborhood_idx(recv_ids(i), j) + if (j > 0) then + ! add forces and subtract recv_snap prevent double-counting + forces(j,:) = forces(j,:) + recv_ft(1:3,i) - recv_forces_snap(j,:) + torques(j,:) = torques(j,:) + recv_ft(4:6,i) - recv_torques_snap(j,:) + recv_forces_snap(j,:) = recv_ft(1:3,i) + recv_torques_snap(j,:) = recv_ft(4:6,i) + end if + end do + $:END_GPU_PARALLEL_LOOP() + end if + tag = tag + 2 + end do + end if + #:endfor + + ! Send final sums back to neighbors in -X direction + #:for X, ID in [('x', 1), ('y', 2), ('z', 3)] + if (num_dims >= ${ID}$) then + send_neighbor = merge(bc_${X}$%beg, MPI_PROC_NULL, bc_${X}$%beg >= 0) + recv_neighbor = merge(bc_${X}$%end, MPI_PROC_NULL, bc_${X}$%end >= 0) + + do k = 1, 2*ib_neighborhood_radius + pack_pos = 0 + $:GPU_PARALLEL_LOOP(private='[i]', copyin='[forces, torques]') + do i = 1, num_ibs + send_ids(i) = patch_ib(i)%gbl_patch_id + send_ft(1:3,i) = forces(i,:) + send_ft(4:6,i) = torques(i,:) + end do + $:END_GPU_PARALLEL_LOOP() + $:GPU_UPDATE(host='[send_ids, send_ft]') + call MPI_PACK(num_ibs, 1, MPI_INTEGER, ib_force_send_buf, buf_size, pack_pos, MPI_COMM_WORLD, ierr) + call MPI_PACK(send_ids, num_ibs, MPI_INTEGER, ib_force_send_buf, buf_size, pack_pos, MPI_COMM_WORLD, ierr) + call MPI_PACK(send_ft, 6*num_ibs, mpi_p, ib_force_send_buf, buf_size, pack_pos, MPI_COMM_WORLD, ierr) + call MPI_SENDRECV(ib_force_send_buf, pack_pos, MPI_PACKED, send_neighbor, tag, ib_force_recv_buf, buf_size, & + & MPI_PACKED, recv_neighbor, tag, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) + if (recv_neighbor /= MPI_PROC_NULL) then + unpack_pos = 0 + call MPI_UNPACK(ib_force_recv_buf, buf_size, unpack_pos, recv_count, 1, MPI_INTEGER, MPI_COMM_WORLD, ierr) + call MPI_UNPACK(ib_force_recv_buf, buf_size, unpack_pos, recv_ids, recv_count, MPI_INTEGER, & + & MPI_COMM_WORLD, ierr) + call MPI_UNPACK(ib_force_recv_buf, buf_size, unpack_pos, recv_ft, 6*recv_count, mpi_p, MPI_COMM_WORLD, ierr) + $:GPU_PARALLEL_LOOP(private='[i, j]', copyin='[recv_ft, recv_ids]', copy='[forces, torques]') + do i = 1, recv_count + call s_get_neighborhood_idx(recv_ids(i), j) + if (j > 0) then + forces(j,:) = recv_ft(1:3,i) + torques(j,:) = recv_ft(4:6,i) + end if + end do + $:END_GPU_PARALLEL_LOOP() + end if + tag = tag + 2 + end do + end if + #:endfor +#endif + + end subroutine s_communicate_ib_forces + + subroutine s_handoff_ib_ownership() + + integer :: i, j, k, output_idx, local_output_idx + integer :: old_num_local_ibs + integer :: new_count, recv_count + integer :: pack_pos, unpack_pos, buf_size, patch_bytes + integer :: send_neighbor, recv_neighbor, ierr + integer :: dx, dy, dz, tag, nbr_idx, nreqs + real(wp), dimension(3) :: centroid + logical :: is_new + type(ib_patch_parameters) :: tmp_patch + integer, dimension(num_local_ibs_max) :: local_ib_idx_old + ! 26 neighbors max in 3D (8 in 2D); each gets its own recv buffer + integer, parameter :: max_nbrs = 26 + character(len=1), allocatable :: send_buf(:), recv_bufs(:,:) + integer, dimension(2*max_nbrs) :: requests + integer, dimension(max_nbrs) :: recv_neighbor_list + +#ifdef MFC_MPI + if (num_procs > 1) then + ! save a copy of the local IB's global indices to cross-reference for later. + local_ib_idx_old = 0 + old_num_local_ibs = num_local_ibs + do i = 1, num_local_ibs + local_ib_idx_old(i) = patch_ib(local_ib_patch_ids(i))%gbl_patch_id + end do + + ! delete any particles that no longer need to be tracked and coalesce the array + output_idx = 0 + local_output_idx = 0 + do i = 1, num_ibs + centroid = [patch_ib(i)%x_centroid, patch_ib(i)%y_centroid, 0._wp] + if (num_dims == 3) centroid(3) = patch_ib(i)%z_centroid + + ! delete if not in neighborhood + if (f_neighborhood_ranks_own_location(centroid)) then + output_idx = output_idx + 1 + if (i /= output_idx) then + patch_ib(output_idx) = patch_ib(i) + end if + + ! check if in local domain + if (f_local_rank_owns_location(centroid)) then + local_output_idx = local_output_idx + 1 + local_ib_patch_ids(local_output_idx) = output_idx + end if + end if + end do + num_ibs = output_idx + num_local_ibs = local_output_idx + $:GPU_UPDATE(device='[patch_ib]') + call s_update_ib_lookup() + + ! Broadcast newly-owned patches to all neighborhood neighbors + patch_bytes = storage_size(tmp_patch)/8 + buf_size = storage_size(0)/8 + patch_bytes*num_local_ibs_max + allocate (send_buf(buf_size), recv_bufs(buf_size, max_nbrs)) + + ! Write placeholder count at position 0 + pack_pos = 0 + call MPI_PACK(0, 1, MPI_INTEGER, send_buf, buf_size, pack_pos, MPI_COMM_WORLD, ierr) + + ! pack new patches and count them + new_count = 0 + do i = 1, num_local_ibs + k = local_ib_patch_ids(i) + is_new = .true. + do j = 1, old_num_local_ibs + if (patch_ib(k)%gbl_patch_id == local_ib_idx_old(j)) then + is_new = .false. + exit + end if + end do + if (is_new) then + call MPI_PACK(patch_ib(k), patch_bytes, MPI_BYTE, send_buf, buf_size, pack_pos, MPI_COMM_WORLD, ierr) + new_count = new_count + 1 + end if + end do + + ! Overwrite the placeholder with the real count + pack_pos = 0 + call MPI_PACK(new_count, 1, MPI_INTEGER, send_buf, buf_size, pack_pos, MPI_COMM_WORLD, ierr) + pack_pos = storage_size(0)/8 + new_count*patch_bytes + + ! Post all receives first, then sends + nreqs = 0 + nbr_idx = 0 + do dz = merge(-1, 0, num_dims == 3), merge(1, 0, num_dims == 3) + do dy = -1, 1 + do dx = -1, 1 + if (dx == 0 .and. dy == 0 .and. dz == 0) cycle + nbr_idx = nbr_idx + 1 + tag = 200 + (dx + 1)*9 + (dy + 1)*3 + (dz + 1) + recv_neighbor = ib_neighbor_ranks(-dx, -dy, -dz) + recv_neighbor_list(nbr_idx) = MPI_PROC_NULL + if (recv_neighbor < 0) cycle + recv_neighbor_list(nbr_idx) = recv_neighbor + nreqs = nreqs + 1 + call MPI_IRECV(recv_bufs(:,nbr_idx), buf_size, MPI_PACKED, recv_neighbor, tag, MPI_COMM_WORLD, & + & requests(nreqs), ierr) + end do + end do + end do + + do dz = merge(-1, 0, num_dims == 3), merge(1, 0, num_dims == 3) + do dy = -1, 1 + do dx = -1, 1 + if (dx == 0 .and. dy == 0 .and. dz == 0) cycle + tag = 200 + (dx + 1)*9 + (dy + 1)*3 + (dz + 1) + send_neighbor = ib_neighbor_ranks(dx, dy, dz) + if (send_neighbor < 0) cycle + nreqs = nreqs + 1 + call MPI_ISEND(send_buf, pack_pos, MPI_PACKED, send_neighbor, tag, MPI_COMM_WORLD, requests(nreqs), ierr) + end do + end do + end do + + call MPI_WAITALL(nreqs, requests, MPI_STATUSES_IGNORE, ierr) + + ! Unpack all received buffers + do nbr_idx = 1, merge(26, 8, num_dims == 3) + if (recv_neighbor_list(nbr_idx) == MPI_PROC_NULL) cycle + unpack_pos = 0 + call MPI_UNPACK(recv_bufs(:,nbr_idx), buf_size, unpack_pos, recv_count, 1, MPI_INTEGER, MPI_COMM_WORLD, ierr) + do i = 1, recv_count + call MPI_UNPACK(recv_bufs(:,nbr_idx), buf_size, unpack_pos, tmp_patch, patch_bytes, MPI_BYTE, MPI_COMM_WORLD, & + & ierr) + call s_get_neighborhood_idx(tmp_patch%gbl_patch_id, j) + if (j < 0) then + num_ibs = num_ibs + 1 + @:ASSERT(num_ibs <= size(patch_ib), 'patch_ib overflow in neighborhood handoff') + patch_ib(num_ibs) = tmp_patch + end if + end do + end do + + deallocate (send_buf, recv_bufs) + $:GPU_UPDATE(device='[patch_ib]') + call s_update_ib_lookup() + end if +#endif + + end subroutine s_handoff_ib_ownership + + subroutine s_get_neighborhood_idx(gbl_idx, neighborhood_idx) + + $:GPU_ROUTINE(parallelism='[seq]') + + integer, intent(in) :: gbl_idx + integer, intent(out) :: neighborhood_idx + integer :: i + + neighborhood_idx = ib_gbl_idx_lookup(gbl_idx) + + end subroutine s_get_neighborhood_idx + + subroutine s_update_ib_lookup() + + integer :: i + + ib_gbl_idx_lookup = -1 + $:GPU_UPDATE(device='[ib_gbl_idx_lookup]') + + $:GPU_PARALLEL_LOOP(private='[i]') + do i = 1, num_ibs + ib_gbl_idx_lookup(patch_ib(i)%gbl_patch_id) = i + end do + $:END_GPU_PARALLEL_LOOP() + + $:GPU_UPDATE(host='[ib_gbl_idx_lookup]') + + end subroutine s_update_ib_lookup + + !> Finalize the IBM module + impure subroutine s_finalize_ibm_module() + + @:DEALLOCATE(ib_markers%sf) + @:DEALLOCATE(ib_gbl_idx_lookup) + if (allocated(airfoil_grid_u)) then + @:DEALLOCATE(airfoil_grid_u) + @:DEALLOCATE(airfoil_grid_l) + end if + + if (allocated(models)) then + @:DEALLOCATE(models) + end if + + if (collision_model > 0) call s_finalize_collisions_module() + +#ifdef MFC_MPI + if (num_procs > 1) then + @:DEALLOCATE(send_ids, send_ft) + deallocate (recv_forces_snap, recv_torques_snap, recv_ids, recv_ft) + end if +#endif + + end subroutine s_finalize_ibm_module + end module m_ibm diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp index b3c7cb6f49..9a7d7326de 100644 --- a/src/simulation/m_mpi_proxy.fpp +++ b/src/simulation/m_mpi_proxy.fpp @@ -74,9 +74,10 @@ contains & 'wave_speeds', 'avg_state', 'precision', 'bc_x%beg', 'bc_x%end', & & 'bc_y%beg', 'bc_y%end', 'bc_z%beg', 'bc_z%end', 'fd_order', & & 'num_probes', 'num_integrals', 'bubble_model', 'thermal', & - & 'num_source', 'relax_model', 'num_ibs', 'n_start', & + & 'num_source', 'relax_model', 'num_ibs', 'num_particle_beds', 'n_start', & & 'num_bc_patches', 'num_igr_iters', 'num_igr_warm_start_iters', & - & 'adap_dt_max_iters', 'int_comp', 'collision_model' ] + & 'adap_dt_max_iters', 'collision_model', 'ib_neighborhood_radius', & + & 'int_comp' ] call MPI_BCAST(${VAR}$, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) #:endfor @@ -191,6 +192,21 @@ contains #:endfor end do + ! Non-rank-0 processes may have patch_ib sized to num_ib_patches_max_namelist while rank 0 grew it + ! for particle beds. Grow here before receiving the broadcast entries. + if (proc_rank /= 0 .and. num_ibs > size(patch_ib)) then + block + type(ib_patch_parameters), allocatable :: tmp(:) + integer :: n + n = size(patch_ib) + allocate (tmp(n)) + tmp(1:n) = patch_ib(1:n) + deallocate (patch_ib) + allocate (patch_ib(num_ib_patches_max)) + patch_ib(1:n) = tmp + end block + end if + do i = 1, num_ibs #:for VAR in [ 'radius', 'length_x', 'length_y', 'length_z', & & 'x_centroid', 'y_centroid', 'z_centroid', 'c', 'm', 'p', 't', 'theta', 'slip', 'mass', & @@ -206,6 +222,16 @@ contains call MPI_BCAST(patch_ib(i)%model_filepath, len(patch_ib(i)%model_filepath), MPI_CHARACTER, 0, MPI_COMM_WORLD, ierr) end do + do i = 1, num_particle_beds + #:for VAR in ['x_centroid', 'y_centroid', 'z_centroid', 'length_x', 'length_y', 'length_z', & + & 'radius', 'mass', 'min_spacing'] + call MPI_BCAST(particle_bed(i)%${VAR}$, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) + #:endfor + call MPI_BCAST(particle_bed(i)%num_particles, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(particle_bed(i)%moving_ibm, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(particle_bed(i)%seed, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + end do + do j = 1, num_probes_max do i = 1, 3 call MPI_BCAST(acoustic(j)%loc(i), 1, mpi_p, 0, MPI_COMM_WORLD, ierr) diff --git a/src/simulation/m_particle_bed.fpp b/src/simulation/m_particle_bed.fpp new file mode 100644 index 0000000000..d8e383ca07 --- /dev/null +++ b/src/simulation/m_particle_bed.fpp @@ -0,0 +1,197 @@ +!> +!! @file m_particle_bed.fpp +!! @brief Generates particle beds: converts particle_bed specifications into +!! individual sphere/circle patch_ib entries before MPI broadcast. + +!> @brief Generates particle beds by converting particle_bed patch specifications into individual immersed boundary patches before +!! MPI broadcast. +module m_particle_bed + + use m_global_parameters + use m_constants + + implicit none + + private + + public :: s_generate_particle_beds + +contains + + !> Generate all particle beds and append the resulting particles to patch_ib. Called on rank 0 only, before + !! s_mpi_bcast_user_inputs. Uses a spatial hash grid (cell size = min_dist) so each candidate requires only 3^dim distance + !! checks on average instead of O(n). + impure subroutine s_generate_particle_beds() + + integer :: b, ib_idx, geom + integer :: n_placed + integer(8) :: n_attempts, max_attempts + real(wp) :: xmin, xmax, ymin, ymax, zmin, zmax, min_dist + real(wp) :: rx, ry, rz, dist + integer :: seed + logical :: overlaps + real(wp), allocatable :: placed(:,:) + + ! Spatial hash grid + integer :: hash_size, slot + integer :: bx, by, bz, nbx, nby, nbz + integer :: dx_b, dy_b, dz_b, dz_lo, dz_hi, j + integer, allocatable :: hash_head(:), chain_next(:) + + if (num_particle_beds == 0) return + + ! Grow patch_ib from the namelist-sized allocation to the full capacity needed for particle beds + if (size(patch_ib) < num_ib_patches_max) then + block + type(ib_patch_parameters), allocatable :: tmp(:) + integer :: n + n = size(patch_ib) + allocate (tmp(n)) + tmp(1:n) = patch_ib(1:n) + deallocate (patch_ib) + allocate (patch_ib(num_ib_patches_max)) + patch_ib(1:n) = tmp + end block + end if + + do b = 1, num_particle_beds + xmin = particle_bed(b)%x_centroid - 0.5_wp*particle_bed(b)%length_x + xmax = particle_bed(b)%x_centroid + 0.5_wp*particle_bed(b)%length_x + ymin = particle_bed(b)%y_centroid - 0.5_wp*particle_bed(b)%length_y + ymax = particle_bed(b)%y_centroid + 0.5_wp*particle_bed(b)%length_y + zmin = particle_bed(b)%z_centroid - 0.5_wp*particle_bed(b)%length_z + zmax = particle_bed(b)%z_centroid + 0.5_wp*particle_bed(b)%length_z + + min_dist = 2._wp*particle_bed(b)%radius + particle_bed(b)%min_spacing + + if (p == 0) then + geom = 2 ! circle for 2D + dz_lo = 0 + dz_hi = 0 + else + geom = 8 ! sphere for 3D + dz_lo = -1 + dz_hi = 1 + end if + + max_attempts = int(particle_bed(b)%num_particles, 8)*10000_8 + n_placed = 0 + n_attempts = 0 + seed = particle_bed(b)%seed + if (seed == 0) seed = 1 + b*1013904223 + + allocate (placed(3, particle_bed(b)%num_particles)) + + ! Hash table: 4x overprovisioned for ~25% load factor, minimum 16 buckets. chain_next(i) links placed particle i to the + ! previous occupant of its bucket. + hash_size = max(16, 4*particle_bed(b)%num_particles) + allocate (hash_head(hash_size)) + allocate (chain_next(particle_bed(b)%num_particles)) + hash_head = -1 + chain_next = -1 + + do while (n_placed < particle_bed(b)%num_particles .and. n_attempts < max_attempts) + n_attempts = n_attempts + 1 + + rx = xmin + f_xorshift(seed)*(xmax - xmin) + ry = ymin + f_xorshift(seed)*(ymax - ymin) + if (p == 0) then + rz = particle_bed(b)%z_centroid + else + rz = zmin + f_xorshift(seed)*(zmax - zmin) + end if + + bx = int(floor(rx/min_dist)) + by = int(floor(ry/min_dist)) + bz = 0 + if (p /= 0) bz = int(floor(rz/min_dist)) + + ! Check 3x3(x3) neighboring bins - O(1) average via hash lookup + overlaps = .false. + outer: do dx_b = -1, 1 + do dy_b = -1, 1 + do dz_b = dz_lo, dz_hi + nbx = bx + dx_b + nby = by + dy_b + nbz = bz + dz_b + slot = f_bin_hash(nbx, nby, nbz, hash_size) + j = hash_head(slot) + do while (j > 0) + if (p == 0) then + dist = sqrt((rx - placed(1, j))**2 + (ry - placed(2, j))**2) + else + dist = sqrt((rx - placed(1, j))**2 + (ry - placed(2, j))**2 + (rz - placed(3, j))**2) + end if + if (dist < min_dist) then + overlaps = .true. + exit outer + end if + j = chain_next(j) + end do + end do + end do + end do outer + + if (.not. overlaps) then + n_placed = n_placed + 1 + placed(1, n_placed) = rx + placed(2, n_placed) = ry + placed(3, n_placed) = rz + + ! Insert into hash grid as head of bucket chain + slot = f_bin_hash(bx, by, bz, hash_size) + chain_next(n_placed) = hash_head(slot) + hash_head(slot) = n_placed + + num_ibs = num_ibs + 1 + ib_idx = num_ibs + + patch_ib(ib_idx)%gbl_patch_id = ib_idx + patch_ib(ib_idx)%geometry = geom + patch_ib(ib_idx)%x_centroid = rx + patch_ib(ib_idx)%y_centroid = ry + patch_ib(ib_idx)%z_centroid = rz + patch_ib(ib_idx)%radius = particle_bed(b)%radius + patch_ib(ib_idx)%mass = particle_bed(b)%mass + patch_ib(ib_idx)%moving_ibm = particle_bed(b)%moving_ibm + end if + end do + + if (n_placed < particle_bed(b)%num_particles) then + print '("WARNING: particle_bed(",I0,"): placed ",I0," of ",I0," particles after ",I0," attempts")', b, n_placed, & + & particle_bed(b)%num_particles, n_attempts + end if + + deallocate (placed, hash_head, chain_next) + end do + + end subroutine s_generate_particle_beds + + !> Xorshift PRNG. Advances seed in-place and returns a value in [0, 1). + function f_xorshift(seed) result(rval) + + integer, intent(inout) :: seed + real(wp) :: rval + + seed = ieor(seed, ishft(seed, 13)) + seed = ieor(seed, ishft(seed, -17)) + seed = ieor(seed, ishft(seed, 5)) + + rval = abs(real(seed, wp))/real(huge(seed), wp) + + end function f_xorshift + + !> Hash bin coordinates to a 1-indexed slot in [1, hash_size]. Uses large prime multipliers to spread bins across buckets. Hash + !! collisions are benign: the distance check catches false neighbours. + function f_bin_hash(bx, by, bz, hash_size) result(slot) + + integer, intent(in) :: bx, by, bz, hash_size + integer :: slot + integer(8) :: key + + key = ieor(ieor(int(bx, 8)*73856093_8, int(by, 8)*19349663_8), int(bz, 8)*83492791_8) + slot = int(mod(abs(key), int(hash_size, 8))) + 1 + + end function f_bin_hash + +end module m_particle_bed diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 65cc825e40..59d1b39507 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -40,6 +40,9 @@ module m_start_up use m_nvtx use m_ibm + use m_model + use m_particle_bed + use m_collisions use m_compile_specific use m_checker_common use m_checker @@ -91,7 +94,8 @@ contains x_a, y_a, z_a, x_b, y_b, z_b, & x_domain, y_domain, z_domain, & hypoelasticity, & - ib, num_ibs, patch_ib, & + ib, num_ibs, ib_neighborhood_radius, patch_ib, & + num_particle_beds, particle_bed, & collision_model, coefficient_of_restitution, collision_time, & ib_coefficient_of_friction, ib_state_wrt, & fluid_pp, bub_pp, probe_wrt, prim_vars_wrt, & @@ -919,6 +923,8 @@ contains else if (t_step_start > 0) then call s_read_ib_restart_data(t_step_start) end if + call s_instantiate_STL_models() + call s_reduce_ib_patch_array() call s_ibm_setup() if (t_step_start == 0 .or. (cfl_dt .and. n_start == 0)) then call s_write_ib_data_file(0) @@ -999,6 +1005,7 @@ contains call s_assign_default_values_to_user_inputs() call s_read_input_file() call s_check_input_file() + call s_generate_particle_beds() print '(" Simulating a ", A, " ", I0, "x", I0, "x", I0, " case on ", I0, " rank(s) ", A, ".")', & #:if not MFC_CASE_OPTIMIZATION @@ -1014,6 +1021,8 @@ contains #else "on CPUs" #endif + else + allocate (patch_ib(num_ib_patches_max_namelist)) end if call s_mpi_bcast_user_inputs() @@ -1197,4 +1206,283 @@ contains end subroutine s_read_ib_restart_data + !> @brief takes the patch_ib struct array that contains all global IB patches and reduces to only contain patches that are in + !! the local computational domain. + subroutine s_reduce_ib_patch_array() + + type(ib_patch_parameters), allocatable :: patch_ib_gbl(:) + real(wp), dimension(3) :: centroid + integer :: i, j + integer :: num_aware_ibs + logical :: is_in_neighborhood, is_local + + ! do all set up for moving immersed boundaries + + moving_immersed_boundary_flag = .false. + do i = 1, num_ibs + if (patch_ib(i)%moving_ibm /= 0) then + moving_immersed_boundary_flag = .true. + exit + end if + end do + + allocate (patch_ib_gbl(num_ibs)) + patch_ib_gbl(1:num_ibs) = patch_ib(1:num_ibs) + call get_neighbor_bounds() ! make sure the bounds of the neighbors are correctly set up + call s_compute_ib_neighbor_ranks() ! build lookup of all neighbor MPI ranks + + num_gbl_ibs = num_ibs + num_local_ibs = num_ibs + do i = 1, num_local_ibs_max + local_ib_patch_ids(i) = i + end do + + deallocate (patch_ib) + +#ifdef MFC_MPI + if (num_procs == 1) then + ! single-rank: every patch is local; allocate to exact size and copy + allocate (patch_ib(num_gbl_ibs)) + patch_ib(1:num_gbl_ibs) = patch_ib_gbl(1:num_gbl_ibs) + deallocate (patch_ib_gbl) + else + ! multi-rank: carve out the local neighbourhood subset + num_aware_ibs = min(num_local_ibs_max*(2*ib_neighborhood_radius + 1)**num_dims, num_ib_patches_max) + allocate (patch_ib(num_aware_ibs)) + + num_local_ibs = 0 + num_ibs = 0 + do i = 1, num_gbl_ibs + ! catch the edge case where th collision lies just outside the computational domain + is_in_neighborhood = .true. + is_local = .true. + centroid = [patch_ib_gbl(i)%x_centroid, patch_ib_gbl(i)%y_centroid, 0._wp] + if (num_dims == 3) centroid(3) = patch_ib_gbl(i)%z_centroid + + if (f_neighborhood_ranks_own_location(centroid)) then + num_ibs = num_ibs + 1 + patch_ib(num_ibs) = patch_ib_gbl(i) + patch_ib(num_ibs)%gbl_patch_id = i + if (f_local_rank_owns_location(centroid)) then + num_local_ibs = num_local_ibs + 1 + local_ib_patch_ids(num_local_ibs) = num_ibs + end if + end if + end do + + deallocate (patch_ib_gbl) + end if +#else + ! no-MPI: every patch is local; allocate to exact size and copy + allocate (patch_ib(num_gbl_ibs)) + patch_ib(1:num_gbl_ibs) = patch_ib_gbl(1:num_gbl_ibs) + deallocate (patch_ib_gbl) +#endif + + @:ALLOCATE(ib_gbl_idx_lookup(1:num_gbl_ibs)) + + end subroutine s_reduce_ib_patch_array + + !> Build ib_neighbor_ranks(-1:1,-1:1,-1:1): MPI ranks of all neighbor domains. Uses two rounds of MPI_SENDRECV cascades - face + !! neighbors are known from bc_*, edge neighbors are obtained in round 1, and (3D) corner neighbors in round 2. + subroutine s_compute_ib_neighbor_ranks() + + integer :: ax, k, nbr_idx, nreqs, sx, sy, sz, dx, dy, dz + integer, allocatable :: send_table(:,:,:), recv_tables(:,:,:,:) + integer, dimension(52) :: requests + +#ifdef MFC_MPI + integer :: ierr + integer, dimension(4) :: buf4, rbuf4 + integer, dimension(2) :: buf2, rbuf2 + + ax = ib_neighborhood_radius + + if (allocated(ib_neighbor_ranks)) deallocate (ib_neighbor_ranks) + allocate (ib_neighbor_ranks(-ax:ax,-ax:ax,-ax:ax)) + ib_neighbor_ranks = MPI_PROC_NULL + ib_neighbor_ranks(0, 0, 0) = proc_rank + + ! Fill radius-1 entries: face neighbors are known from domain decomposition + ib_neighbor_ranks(-1, 0, 0) = bc_x%beg + ib_neighbor_ranks(+1, 0, 0) = bc_x%end + if (num_dims >= 2) then + ib_neighbor_ranks(0, -1, 0) = bc_y%beg + ib_neighbor_ranks(0, +1, 0) = bc_y%end + end if + if (num_dims == 3) then + ib_neighbor_ranks(0, 0, -1) = bc_z%beg + ib_neighbor_ranks(0, 0, +1) = bc_z%end + end if + + if (num_dims >= 2) then + ! Round 1a: exchange y/z face ranks with +/-x face neighbors -> xy and xz edge ranks + buf4 = [bc_y%beg, bc_y%end, bc_z%beg, bc_z%end] + + ! Send to -x, receive from +x -> edges (+1,+/-1,0) and (+1,0,+/-1) + call MPI_SENDRECV(buf4, 4, MPI_INTEGER, merge(bc_x%beg, MPI_PROC_NULL, bc_x%beg >= 0), 310, rbuf4, 4, MPI_INTEGER, & + & merge(bc_x%end, MPI_PROC_NULL, bc_x%end >= 0), 310, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) + if (bc_x%end >= 0) then + ib_neighbor_ranks(+1, -1, 0) = rbuf4(1) + ib_neighbor_ranks(+1, +1, 0) = rbuf4(2) + ib_neighbor_ranks(+1, 0, -1) = rbuf4(3) + ib_neighbor_ranks(+1, 0, +1) = rbuf4(4) + end if + + call MPI_SENDRECV(buf4, 4, MPI_INTEGER, merge(bc_x%end, MPI_PROC_NULL, bc_x%end >= 0), 311, rbuf4, 4, MPI_INTEGER, & + & merge(bc_x%beg, MPI_PROC_NULL, bc_x%beg >= 0), 311, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) + if (bc_x%beg >= 0) then + ib_neighbor_ranks(-1, -1, 0) = rbuf4(1) + ib_neighbor_ranks(-1, +1, 0) = rbuf4(2) + ib_neighbor_ranks(-1, 0, -1) = rbuf4(3) + ib_neighbor_ranks(-1, 0, +1) = rbuf4(4) + end if + end if + + if (num_dims == 3) then + ! Round 1b: exchange z face ranks with +/-y face neighbors -> yz edge ranks + buf2 = [bc_z%beg, bc_z%end] + + call MPI_SENDRECV(buf2, 2, MPI_INTEGER, merge(bc_y%beg, MPI_PROC_NULL, bc_y%beg >= 0), 312, rbuf2, 2, MPI_INTEGER, & + & merge(bc_y%end, MPI_PROC_NULL, bc_y%end >= 0), 312, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) + if (bc_y%end >= 0) then + ib_neighbor_ranks(0, +1, -1) = rbuf2(1) + ib_neighbor_ranks(0, +1, +1) = rbuf2(2) + end if + + call MPI_SENDRECV(buf2, 2, MPI_INTEGER, merge(bc_y%end, MPI_PROC_NULL, bc_y%end >= 0), 313, rbuf2, 2, MPI_INTEGER, & + & merge(bc_y%beg, MPI_PROC_NULL, bc_y%beg >= 0), 313, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) + if (bc_y%beg >= 0) then + ib_neighbor_ranks(0, -1, -1) = rbuf2(1) + ib_neighbor_ranks(0, -1, +1) = rbuf2(2) + end if + + ! Round 2: exchange z face ranks with xy-diagonal edge neighbors -> corner ranks. Each of the 4 xy diagonals gives 2 + ! corners (the +/-z variants). Pattern: send buf2 to mirror diagonal, receive from this diagonal -> that edge's z face + ! ranks. + #:for DX, DY, MDX, MDY, TAG in [(1,1,-1,-1,320), (1,-1,-1,1,321), (-1,1,1,-1,322), (-1,-1,1,1,323)] + call MPI_SENDRECV(buf2, 2, MPI_INTEGER, merge(ib_neighbor_ranks(${MDX}$, ${MDY}$, 0), MPI_PROC_NULL, & + & ib_neighbor_ranks(${MDX}$, ${MDY}$, 0) >= 0), ${TAG}$, rbuf2, 2, MPI_INTEGER, & + & merge(ib_neighbor_ranks(${DX}$, ${DY}$, 0), MPI_PROC_NULL, ib_neighbor_ranks(${DX}$, ${DY}$, & + & 0) >= 0), ${TAG}$, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) + if (ib_neighbor_ranks(${DX}$, ${DY}$, 0) >= 0) then + ib_neighbor_ranks(${DX}$, ${DY}$, -1) = rbuf2(1) + ib_neighbor_ranks(${DX}$, ${DY}$, +1) = rbuf2(2) + end if + #:endfor + end if + + ! For radius > 1: extend the table by iterative 26-neighbor full-table exchanges. In each round, every rank broadcasts its + ! current table to all 26 immediate neighbors. Their entry at offset (dx,dy,dz) from them = our entry at + ! (dx+sx,dy+sy,dz+sz). One extension round fills the entire next shell, so ax-1 rounds suffice. + if (ax > 1) then + allocate (send_table(-ax:ax,-ax:ax,-ax:ax)) + allocate (recv_tables(-ax:ax,-ax:ax,-ax:ax,1:26)) + + do k = 2, ax + send_table = ib_neighbor_ranks + + nreqs = 0 + nbr_idx = 0 + do sz = -1, 1 + do sy = -1, 1 + do sx = -1, 1 + if (sx == 0 .and. sy == 0 .and. sz == 0) cycle + nbr_idx = nbr_idx + 1 + if (ib_neighbor_ranks(sx, sy, sz) < 0) cycle + nreqs = nreqs + 1 + call MPI_IRECV(recv_tables(:,:,:,nbr_idx), (2*ax + 1)**3, MPI_INTEGER, ib_neighbor_ranks(sx, sy, sz), & + & 400, MPI_COMM_WORLD, requests(nreqs), ierr) + end do + end do + end do + + do sz = -1, 1 + do sy = -1, 1 + do sx = -1, 1 + if (sx == 0 .and. sy == 0 .and. sz == 0) cycle + if (ib_neighbor_ranks(sx, sy, sz) < 0) cycle + nreqs = nreqs + 1 + call MPI_ISEND(send_table, (2*ax + 1)**3, MPI_INTEGER, ib_neighbor_ranks(sx, sy, sz), 400, & + & MPI_COMM_WORLD, requests(nreqs), ierr) + end do + end do + end do + + call MPI_WAITALL(nreqs, requests, MPI_STATUSES_IGNORE, ierr) + + nbr_idx = 0 + do sz = -1, 1 + do sy = -1, 1 + do sx = -1, 1 + if (sx == 0 .and. sy == 0 .and. sz == 0) cycle + nbr_idx = nbr_idx + 1 + if (ib_neighbor_ranks(sx, sy, sz) < 0) cycle + do dz = -ax, ax + do dy = -ax, ax + do dx = -ax, ax + if (recv_tables(dx, dy, dz, nbr_idx) == MPI_PROC_NULL) cycle + if (dx + sx < -ax .or. dx + sx > ax) cycle + if (dy + sy < -ax .or. dy + sy > ax) cycle + if (dz + sz < -ax .or. dz + sz > ax) cycle + if (ib_neighbor_ranks(dx + sx, dy + sy, dz + sz) /= MPI_PROC_NULL) cycle + ib_neighbor_ranks(dx + sx, dy + sy, dz + sz) = recv_tables(dx, dy, dz, nbr_idx) + end do + end do + end do + end do + end do + end do + end do + + deallocate (send_table, recv_tables) + end if +#endif + + end subroutine s_compute_ib_neighbor_ranks + + subroutine get_neighbor_bounds() + + real(wp) :: beg_val, end_val, recv_val + integer :: k, send_neighbor, recv_neighbor, ierr + + ! Default: unbounded in all directions (covers single-rank and no-MPI cases) + + neighbor_domain_x%beg = -huge(0._wp) + neighbor_domain_x%end = huge(0._wp) + neighbor_domain_y%beg = -huge(0._wp) + neighbor_domain_y%end = huge(0._wp) + neighbor_domain_z%beg = -huge(0._wp) + neighbor_domain_z%end = huge(0._wp) + +#ifdef MFC_MPI + ! For each direction, propagate the left/right boundary edges outward ib_neighborhood_radius hops. After k rounds: beg_val = + ! left edge of the rank k hops to the left; end_val = right edge of the rank k hops to the right. + #:for X, ID, TAG, DIM in [('x', 1, 100, 'm'), ('y', 2, 102, 'n'), ('z', 3, 104, 'p')] + if (num_dims >= ${ID}$) then + beg_val = ${X}$_cb(-1) + end_val = ${X}$_cb(${DIM}$) + do k = 1, ib_neighborhood_radius + send_neighbor = merge(bc_${X}$%end, MPI_PROC_NULL, bc_${X}$%end >= 0) + recv_neighbor = merge(bc_${X}$%beg, MPI_PROC_NULL, bc_${X}$%beg >= 0) + recv_val = -huge(0._wp) + call MPI_SENDRECV(beg_val, 1, mpi_p, send_neighbor, ${TAG}$, recv_val, 1, mpi_p, recv_neighbor, ${TAG}$, & + & MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) + beg_val = recv_val + + send_neighbor = merge(bc_${X}$%beg, MPI_PROC_NULL, bc_${X}$%beg >= 0) + recv_neighbor = merge(bc_${X}$%end, MPI_PROC_NULL, bc_${X}$%end >= 0) + recv_val = huge(0._wp) + call MPI_SENDRECV(end_val, 1, mpi_p, send_neighbor, ${TAG + 1}$, recv_val, 1, mpi_p, recv_neighbor, & + & ${TAG + 1}$, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) + end_val = recv_val + end do + neighbor_domain_${X}$%beg = beg_val + neighbor_domain_${X}$%end = end_val + end if + #:endfor +#endif + + end subroutine get_neighbor_bounds + end module m_start_up diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index 130053fce5..66dd1848f8 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -552,7 +552,6 @@ contains ! check if any IBMS are moving, and if so, update the markers, ghost points, levelsets, and levelset norms if (moving_immersed_boundary_flag) then call s_propagate_immersed_boundaries(s) - ! compute ib forces for fixed immersed boundaries if requested for output end if ! update the ghost fluid properties point values based on IB state @@ -564,10 +563,11 @@ contains end if end do - ! if (ib) then - if (moving_immersed_boundary_flag) call s_wrap_periodic_ibs() - if (ib_state_wrt .and. (.not. moving_immersed_boundary_flag)) then + if (moving_immersed_boundary_flag) then + call s_wrap_periodic_ibs() ! wraps the positions of IBs to the local proc + call s_handoff_ib_ownership() ! recomputes which ranks own which IBs and communicate to neighbors + else if (ib_state_wrt) then call s_compute_ib_forces(q_prim_vf, fluid_pp) end if end if @@ -708,11 +708,11 @@ contains integer, intent(in) :: s integer :: i - logical :: forces_computed + integer :: gbl_id ! used for analytic ib patch motion call nvtxStartRange("PROPAGATE-IMMERSED-BOUNDARIES") - forces_computed = .false. + if (moving_immersed_boundary_flag) call s_compute_ib_forces(q_prim_vf, fluid_pp) do i = 1, num_ibs if (s == 1) then @@ -726,11 +726,6 @@ contains ! Compute forces BEFORE the RK velocity blend so the device copy of patch_ib%vel matches the host (pre-blend) when ! velocity-dependent collision damping forces are evaluated on the GPU. - if (patch_ib(i)%moving_ibm == 2 .and. .not. forces_computed) then - call s_compute_ib_forces(q_prim_vf, fluid_pp) - forces_computed = .true. - end if - if (patch_ib(i)%moving_ibm > 0) then patch_ib(i)%vel = (rk_coef(s, 1)*patch_ib(i)%step_vel + rk_coef(s, 2)*patch_ib(i)%vel)/rk_coef(s, 4) patch_ib(i)%angular_vel = (rk_coef(s, 1)*patch_ib(i)%step_angular_vel + rk_coef(s, & @@ -748,7 +743,6 @@ contains & 3)*dt*patch_ib(i)%torque/rk_coef(s, 4)) ! add the torque to the angular momentum call s_compute_moment_of_inertia(i, patch_ib(i)%angular_vel) ! update the moment of inertia to be based on the direction of the angular momentum - ! convert back to angular velocity with the new moment of inertia patch_ib(i)%angular_vel = patch_ib(i)%angular_vel/patch_ib(i)%moment end if diff --git a/tests/BA186DDF/golden-metadata.txt b/tests/BA186DDF/golden-metadata.txt new file mode 100644 index 0000000000..6888112175 --- /dev/null +++ b/tests/BA186DDF/golden-metadata.txt @@ -0,0 +1,190 @@ +This file was created on 2026-05-13 07:35:59.310410. + +mfc.sh: + + Invocation: test --generate --only BA186DDF -j 8 + Lock: mpi=Yes & gpu=No & debug=No & reldebug=No & gcov=No & unified=No & single=No & mixed=No & fastmath=No + Git: 9c2a820a102dc19231f6782d1e4141dc2c482f72 on HEAD (clean) + +pre_process: + + CMake Configuration: + + CMake v3.25.2 on k004-004.hpcfund + + C : GNU v12.2.0 (/opt/ohpc/pub/compiler/gcc/12.2.0/bin/cc) + Fortran : GNU v12.2.0 (/opt/ohpc/pub/compiler/gcc/12.2.0/bin/gfortran) + + PRE_PROCESS : ON + SIMULATION : OFF + POST_PROCESS : OFF + SYSCHECK : OFF + DOCUMENTATION : OFF + ALL : OFF + + MPI : ON + OpenACC : OFF + OpenMP : OFF + + Fypp : /work1/spencerbryngelson/sbryngelson/local-aware-ibm/build/venv/bin/fypp + Doxygen : + + Build Type : Release + + Configuration Environment: + + CC : /opt/ohpc/pub/compiler/gcc/12.2.0/bin/cc + CXX : /opt/ohpc/pub/compiler/gcc/12.2.0/bin/c++ + FC : /opt/ohpc/pub/compiler/gcc/12.2.0/bin/gfortran + OMPI_CC : + OMPI_CXX : + OMPI_FC : + +post_process: + + CMake Configuration: + + CMake v3.25.2 on k004-004.hpcfund + + C : GNU v12.2.0 (/opt/ohpc/pub/compiler/gcc/12.2.0/bin/cc) + Fortran : GNU v12.2.0 (/opt/ohpc/pub/compiler/gcc/12.2.0/bin/gfortran) + + PRE_PROCESS : OFF + SIMULATION : OFF + POST_PROCESS : ON + SYSCHECK : OFF + DOCUMENTATION : OFF + ALL : OFF + + MPI : ON + OpenACC : OFF + OpenMP : OFF + + Fypp : /work1/spencerbryngelson/sbryngelson/local-aware-ibm/build/venv/bin/fypp + Doxygen : + + Build Type : Release + + Configuration Environment: + + CC : /opt/ohpc/pub/compiler/gcc/12.2.0/bin/cc + CXX : /opt/ohpc/pub/compiler/gcc/12.2.0/bin/c++ + FC : /opt/ohpc/pub/compiler/gcc/12.2.0/bin/gfortran + OMPI_CC : + OMPI_CXX : + OMPI_FC : + +simulation: + + CMake Configuration: + + CMake v3.25.2 on k004-004.hpcfund + + C : GNU v12.2.0 (/opt/ohpc/pub/compiler/gcc/12.2.0/bin/cc) + Fortran : GNU v12.2.0 (/opt/ohpc/pub/compiler/gcc/12.2.0/bin/gfortran) + + PRE_PROCESS : OFF + SIMULATION : ON + POST_PROCESS : OFF + SYSCHECK : OFF + DOCUMENTATION : OFF + ALL : OFF + + MPI : ON + OpenACC : OFF + OpenMP : OFF + + Fypp : /work1/spencerbryngelson/sbryngelson/local-aware-ibm/build/venv/bin/fypp + Doxygen : + + Build Type : Release + + Configuration Environment: + + CC : /opt/ohpc/pub/compiler/gcc/12.2.0/bin/cc + CXX : /opt/ohpc/pub/compiler/gcc/12.2.0/bin/c++ + FC : /opt/ohpc/pub/compiler/gcc/12.2.0/bin/gfortran + OMPI_CC : + OMPI_CXX : + OMPI_FC : + +syscheck: + + CMake Configuration: + + CMake v3.25.2 on k004-004.hpcfund + + C : GNU v12.2.0 (/opt/ohpc/pub/compiler/gcc/12.2.0/bin/cc) + Fortran : GNU v12.2.0 (/opt/ohpc/pub/compiler/gcc/12.2.0/bin/gfortran) + + PRE_PROCESS : OFF + SIMULATION : OFF + POST_PROCESS : OFF + SYSCHECK : ON + DOCUMENTATION : OFF + ALL : OFF + + MPI : ON + OpenACC : OFF + OpenMP : OFF + + Fypp : /work1/spencerbryngelson/sbryngelson/local-aware-ibm/build/venv/bin/fypp + Doxygen : + + Build Type : Release + + Configuration Environment: + + CC : /opt/ohpc/pub/compiler/gcc/12.2.0/bin/cc + CXX : /opt/ohpc/pub/compiler/gcc/12.2.0/bin/c++ + FC : /opt/ohpc/pub/compiler/gcc/12.2.0/bin/gfortran + OMPI_CC : + OMPI_CXX : + OMPI_FC : + +CPU: + + CPU Info: + From lscpu + Architecture: x86_64 + CPU op-mode(s): 32-bit, 64-bit + Address sizes: 48 bits physical, 48 bits virtual + Byte Order: Little Endian + CPU(s): 128 + On-line CPU(s) list: 0-127 + Vendor ID: AuthenticAMD + Model name: AMD EPYC 7763 64-Core Processor + CPU family: 25 + Model: 1 + Thread(s) per core: 1 + Core(s) per socket: 64 + Socket(s): 2 + Stepping: 1 + Frequency boost: enabled + CPU max MHz: 3529.0520 + CPU min MHz: 1500.0000 + BogoMIPS: 4891.26 + Flags: fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush mmx fxsr sse sse2 ht syscall nx mmxext fxsr_opt pdpe1gb rdtscp lm constant_tsc rep_good nopl nonstop_tsc cpuid extd_apicid aperfmperf rapl pni pclmulqdq monitor ssse3 fma cx16 pcid sse4_1 sse4_2 movbe popcnt aes xsave avx f16c rdrand lahf_lm cmp_legacy svm extapic cr8_legacy abm sse4a misalignsse 3dnowprefetch osvw ibs skinit wdt tce topoext perfctr_core perfctr_nb bpext perfctr_llc mwaitx cpb cat_l3 cdp_l3 hw_pstate ssbd mba ibrs ibpb stibp vmmcall fsgsbase bmi1 avx2 smep bmi2 erms invpcid cqm rdt_a rdseed adx smap clflushopt clwb sha_ni xsaveopt xsavec xgetbv1 xsaves cqm_llc cqm_occup_llc cqm_mbm_total cqm_mbm_local clzero irperf xsaveerptr rdpru wbnoinvd amd_ppin brs arat npt lbrv svm_lock nrip_save tsc_scale vmcb_clean flushbyasid decodeassists pausefilter pfthreshold v_vmsave_vmload vgif v_spec_ctrl umip pku ospke vaes vpclmulqdq rdpid overflow_recov succor smca fsrm debug_swap + Virtualization: AMD-V + L1d cache: 4 MiB (128 instances) + L1i cache: 4 MiB (128 instances) + L2 cache: 64 MiB (128 instances) + L3 cache: 512 MiB (16 instances) + NUMA node(s): 2 + NUMA node0 CPU(s): 0-63 + NUMA node1 CPU(s): 64-127 + Vulnerability Gather data sampling: Not affected + Vulnerability Itlb multihit: Not affected + Vulnerability L1tf: Not affected + Vulnerability Mds: Not affected + Vulnerability Meltdown: Not affected + Vulnerability Mmio stale data: Not affected + Vulnerability Reg file data sampling: Not affected + Vulnerability Retbleed: Not affected + Vulnerability Spec rstack overflow: Mitigation; Safe RET + Vulnerability Spec store bypass: Mitigation; Speculative Store Bypass disabled via prctl + Vulnerability Spectre v1: Mitigation; usercopy/swapgs barriers and __user pointer sanitization + Vulnerability Spectre v2: Mitigation; Retpolines; IBPB conditional; IBRS_FW; STIBP disabled; RSB filling; PBRSB-eIBRS Not affected; BHI Not affected + Vulnerability Srbds: Not affected + Vulnerability Tsx async abort: Not affected + diff --git a/tests/BA186DDF/golden.txt b/tests/BA186DDF/golden.txt new file mode 100644 index 0000000000..e40f7bf009 --- /dev/null +++ b/tests/BA186DDF/golden.txt @@ -0,0 +1,10 @@ +D/cons.1.00.000000.dat 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 2.6069 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 +D/cons.1.00.000050.dat 2.60690283495417 2.60690283495417 2.60690283495417 2.60690283495417 2.60690283495417 2.60690283495417 2.60690283495417 2.60690283495417 2.60690283495417 2.60690283495417 2.60690283495417 2.60690283495417 2.60690283495417 2.60690283495417 2.60690283495417 2.60690283495417 2.60690283495417 2.60690283495417 2.60690283495417 2.60690283495417 2.60690283495417 2.60690283495417 2.60690283495417 2.60690283495417 2.60690283495417 2.60690283495417 2.60704515690204 2.60704515690204 2.60704515690204 2.60704515690204 2.60704515690204 2.60704515690204 2.60704515690204 2.60704515690204 2.60704515690204 2.60704515690204 2.60704515690204 2.60704515690204 2.60704515690204 2.60704515690204 2.60704515690204 2.60704515690204 2.60704515690204 2.60704515690204 2.60704515690204 2.60704515690204 2.60704515690204 2.60704515690204 2.60704515690204 2.60704515690204 2.60704515690204 2.60704515690204 2.60839351714553 2.60839351714553 2.60839351714553 2.60839351714553 2.60839351714553 2.60839351714553 2.60839351714553 2.60839351714553 2.60839351714553 2.60839351714553 2.60839351714553 2.60839351714553 2.60839351714553 2.60839351714553 2.60839351714553 2.60839351714553 2.60839351714553 2.60839351714553 2.60839351714553 2.60839351714553 2.60839351714553 2.60839351714553 2.60839351714553 2.60839351714553 2.60839351714553 2.60839351714553 1.81733171782674 1.81733171782674 1.81733171782674 1.81733171782674 1.81733171782674 1.81733171782674 1.81733171782674 1.81733171782674 1.81733171782674 1.81733171782674 1.81733171782674 1.81733171782674 1.81733171782674 1.81733171782674 1.81733171782674 1.81733171782674 1.81733171782667 1.81733171782664 1.81733171782674 1.81733171782674 1.81733171782674 1.81733171782674 1.81733171782674 1.81733171782674 1.81733171782674 1.81733171782674 1.42165257910932 1.42165257910932 1.42165257910932 1.42165257910932 1.42165257910932 1.42165257910932 1.42165257910932 1.42165257910932 1.42165257910932 1.42165257910932 1.42165257910932 1.42165257910932 1.42165257910932 1.42165257910932 1.42165257910932 1.42165257910932 1.4216525791091 1.42165257910899 1.42165257910932 1.42165257910932 1.42165257910932 1.42165257910932 1.42165257910932 1.42165257910932 1.42165257910932 1.42165257910932 1.40060529636878 1.40060529636878 1.40060529636878 1.40060529636878 1.40060529636878 1.40060529636878 1.40060529636878 1.40060529636878 1.40060529636878 1.40060529636878 1.40060529636878 1.40060529636878 1.40060529636878 1.40060529636878 1.40060529636878 1.40060529636881 1.40060529637106 1.40060529637204 1.40060529636887 1.40060529636878 1.40060529636878 1.40060529636878 1.40060529636878 1.40060529636878 1.40060529636878 1.40060529636878 1.4000125435324 1.4000125435324 1.4000125435324 1.4000125435324 1.4000125435324 1.4000125435324 1.4000125435324 1.4000125435324 1.4000125435324 1.4000125435324 1.4000125435324 1.4000125435324 1.4000125435324 1.4000125435324 1.4000125435324 1.40001254353241 1.40001254353394 1.40001254353476 1.40001254353244 1.4000125435324 1.4000125435324 1.4000125435324 1.4000125435324 1.4000125435324 1.4000125435324 1.4000125435324 1.40000021079562 1.40000021079562 1.40000021079562 1.40000021079562 1.40000021079562 1.40000021079562 1.40000021079562 1.40000021079562 1.40000021079562 1.40000021079562 1.40000021079562 1.40000021079562 1.40000021079562 1.40000021079562 1.40000021079562 1.40000021079527 1.40000021077065 1.4000002107594 1.40000021079468 1.40000021079561 1.40000021079562 1.40000021079562 1.40000021079562 1.40000021079562 1.40000021079562 1.40000021079562 1.40000000296246 1.40000000296246 1.40000000296246 1.40000000296246 1.40000000296246 1.40000000296246 1.40000000296246 1.40000000296246 1.40000000296246 1.40000000296246 1.40000000296246 1.40000000296246 1.40000000296246 1.40000000296246 1.40000000296246 1.40000000296318 1.40000000300135 1.40000000301236 1.40000000296411 1.40000000296248 1.40000000296246 1.40000000296246 1.40000000296246 1.40000000296246 1.40000000296246 1.40000000296246 1.40000000003006 1.40000000003006 1.40000000003006 1.40000000003006 1.40000000003006 1.40000000003006 1.40000000003006 1.40000000003006 1.40000000003006 1.40000000003006 1.40000000003006 1.40000000003006 1.40000000003006 1.40000000003007 1.40000000003023 1.40000000004039 1.40000000126991 1.40000000183244 1.40000000005878 1.40000000003054 1.40000000003007 1.40000000003006 1.40000000003006 1.40000000003006 1.40000000003006 1.40000000003006 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.40000000000001 1.4 1.40000000000001 1.40000000000003 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 +D/cons.2.00.000000.dat 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 1.81023136 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 +D/cons.2.00.000050.dat 1.81023007117005 1.81023007117005 1.81023007117005 1.81023007117005 1.81023007117005 1.81023007117005 1.81023007117005 1.81023007117005 1.81023007117005 1.81023007117005 1.81023007117005 1.81023007117005 1.81023007117005 1.81023007117005 1.81023007117005 1.81023007117005 1.81023007117005 1.81023007117005 1.81023007117005 1.81023007117005 1.81023007117005 1.81023007117005 1.81023007117005 1.81023007117005 1.81023007117005 1.81023007117005 1.81016539579586 1.81016539579586 1.81016539579586 1.81016539579586 1.81016539579586 1.81016539579586 1.81016539579586 1.81016539579586 1.81016539579586 1.81016539579586 1.81016539579586 1.81016539579586 1.81016539579586 1.81016539579586 1.81016539579586 1.81016539579586 1.81016539579586 1.81016539579586 1.81016539579586 1.81016539579586 1.81016539579586 1.81016539579586 1.81016539579586 1.81016539579586 1.81016539579586 1.81016539579586 1.80972917206017 1.80972917206017 1.80972917206017 1.80972917206017 1.80972917206017 1.80972917206017 1.80972917206017 1.80972917206017 1.80972917206017 1.80972917206017 1.80972917206017 1.80972917206017 1.80972917206017 1.80972917206017 1.80972917206017 1.80972917206017 1.80972917206016 1.80972917206016 1.80972917206017 1.80972917206017 1.80972917206017 1.80972917206017 1.80972917206017 1.80972917206017 1.80972917206017 1.80972917206017 0.63766622233616 0.63766622233616 0.63766622233616 0.63766622233616 0.63766622233616 0.63766622233616 0.63766622233616 0.63766622233616 0.63766622233616 0.63766622233616 0.63766622233616 0.63766622233616 0.63766622233616 0.63766622233616 0.63766622233616 0.63766622233616 0.63766622233624 0.63766622233627 0.63766622233616 0.63766622233616 0.63766622233616 0.63766622233616 0.63766622233616 0.63766622233616 0.63766622233616 0.63766622233616 0.02414334818285 0.02414334818285 0.02414334818285 0.02414334818285 0.02414334818285 0.02414334818285 0.02414334818285 0.02414334818285 0.02414334818285 0.02414334818285 0.02414334818285 0.02414334818285 0.02414334818285 0.02414334818285 0.02414334818285 0.02414334818286 0.02414334818307 0.02414334818318 0.02414334818286 0.02414334818285 0.02414334818285 0.02414334818285 0.02414334818285 0.02414334818285 0.02414334818285 0.02414334818285 0.00060751246729 0.00060751246729 0.00060751246729 0.00060751246729 0.00060751246729 0.00060751246729 0.00060751246729 0.00060751246729 0.00060751246729 0.00060751246729 0.00060751246729 0.00060751246729 0.00060751246729 0.00060751246729 0.00060751246729 0.00060751246727 0.00060751246505 0.00060751246408 0.00060751246723 0.00060751246729 0.00060751246729 0.00060751246729 0.00060751246729 0.00060751246729 0.00060751246729 0.00060751246729 1.254481946e-05 1.254481946e-05 1.254481946e-05 1.254481946e-05 1.254481946e-05 1.254481946e-05 1.254481946e-05 1.254481946e-05 1.254481946e-05 1.254481946e-05 1.254481946e-05 1.254481946e-05 1.254481946e-05 1.254481946e-05 1.254481946e-05 1.254481946e-05 1.254481735e-05 1.254481624e-05 1.254481947e-05 1.254481946e-05 1.254481946e-05 1.254481946e-05 1.254481946e-05 1.254481946e-05 1.254481946e-05 1.254481946e-05 2.1080499e-07 2.1080499e-07 2.1080499e-07 2.1080499e-07 2.1080499e-07 2.1080499e-07 2.1080499e-07 2.1080499e-07 2.1080499e-07 2.1080499e-07 2.1080499e-07 2.1080499e-07 2.1080499e-07 2.1080499e-07 2.1080499e-07 2.1080518e-07 2.1083956e-07 2.1085554e-07 2.1080553e-07 2.10805e-07 2.1080499e-07 2.1080499e-07 2.1080499e-07 2.1080499e-07 2.1080499e-07 2.1080499e-07 2.95347e-09 2.95347e-09 2.95347e-09 2.95347e-09 2.95347e-09 2.95347e-09 2.95347e-09 2.95347e-09 2.95347e-09 2.95347e-09 2.95347e-09 2.95347e-09 2.95347e-09 2.95347e-09 2.95346e-09 2.95271e-09 2.8665e-09 2.83324e-09 2.95134e-09 2.95344e-09 2.95347e-09 2.95347e-09 2.95347e-09 2.95347e-09 2.95347e-09 2.95347e-09 3.007e-11 3.007e-11 3.007e-11 3.007e-11 3.007e-11 3.007e-11 3.007e-11 3.007e-11 3.007e-11 3.007e-11 3.007e-11 3.007e-11 3.007e-11 3.007e-11 3.007e-11 3.064e-11 0.0 0.0 3.166e-11 3.009e-11 3.007e-11 3.007e-11 3.007e-11 3.007e-11 3.007e-11 3.007e-11 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1e-14 0.0 0.0 3e-14 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 0.0 0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 0.0 -0.0 -0.0 0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0 -0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 -0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 0.0 -0.0 0.0 -0.0 -0.0 -0.0 0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 -0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 +D/cons.3.00.000000.dat 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 +D/cons.3.00.000050.dat 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 -0.0 0.0 0.0 -0.0 0.0 0.0 -0.0 0.0 0.0 0.0 -0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 0.0 0.0 -0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 -0.0 0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 -0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 0.0 0.0 -0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 -0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 -0.0 0.0 0.0 0.0 0.0 1e-14 -1e-14 -1e-14 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 -0.0 -2e-14 -8e-14 5e-14 4e-14 0.0 0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -1e-14 -7e-14 4e-14 4e-14 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 2.3e-13 1.03e-12 -6.3e-13 -6.3e-13 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -3.4e-13 -1.32e-12 1.08e-12 5.8e-13 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -1.7e-13 -1.125e-11 0.0 0.0 3.13e-11 4.7e-13 1e-14 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 -0.0 0.0 0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 -0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 -0.0 0.0 0.0 -0.0 0.0 0.0 -0.0 -0.0 0.0 -0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 0.0 -0.0 0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 0.0 -0.0 0.0 -0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 -0.0 0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 -0.0 0.0 -0.0 0.0 0.0 -0.0 0.0 -0.0 0.0 0.0 -0.0 0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 -0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0 0.0 -0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 -0.0 0.0 0.0 -0.0 0.0 -0.0 0.0 -0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 0.0 -0.0 0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 +D/cons.4.00.000000.dat 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 6.774262328192 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 +D/cons.4.00.000050.dat 6.7742701065218 6.7742701065218 6.7742701065218 6.7742701065218 6.7742701065218 6.7742701065218 6.7742701065218 6.7742701065218 6.7742701065218 6.7742701065218 6.7742701065218 6.7742701065218 6.7742701065218 6.7742701065218 6.7742701065218 6.7742701065218 6.7742701065218 6.7742701065218 6.7742701065218 6.7742701065218 6.7742701065218 6.7742701065218 6.7742701065218 6.7742701065218 6.7742701065218 6.7742701065218 6.77466075392754 6.77466075392754 6.77466075392754 6.77466075392754 6.77466075392754 6.77466075392754 6.77466075392754 6.77466075392754 6.77466075392754 6.77466075392754 6.77466075392754 6.77466075392754 6.77466075392754 6.77466075392754 6.77466075392754 6.77466075392754 6.77466075392754 6.77466075392754 6.77466075392754 6.77466075392754 6.77466075392754 6.77466075392754 6.77466075392754 6.77466075392754 6.77466075392754 6.77466075392754 6.77863494620395 6.77863494620395 6.77863494620395 6.77863494620395 6.77863494620395 6.77863494620395 6.77863494620395 6.77863494620395 6.77863494620395 6.77863494620395 6.77863494620395 6.77863494620395 6.77863494620395 6.77863494620395 6.77863494620395 6.77863494620395 6.77863494620395 6.77863494620395 6.77863494620395 6.77863494620395 6.77863494620395 6.77863494620395 6.77863494620395 6.77863494620395 6.77863494620395 6.77863494620395 3.9989398857781 3.9989398857781 3.9989398857781 3.9989398857781 3.9989398857781 3.9989398857781 3.9989398857781 3.9989398857781 3.9989398857781 3.9989398857781 3.9989398857781 3.9989398857781 3.9989398857781 3.9989398857781 3.9989398857781 3.99893988577809 3.99893988577792 3.99893988577784 3.99893988577808 3.9989398857781 3.9989398857781 3.9989398857781 3.9989398857781 3.9989398857781 3.9989398857781 3.9989398857781 2.55743731427511 2.55743731427511 2.55743731427511 2.55743731427511 2.55743731427511 2.55743731427511 2.55743731427511 2.55743731427511 2.55743731427511 2.55743731427511 2.55743731427511 2.55743731427511 2.55743731427511 2.55743731427511 2.55743731427511 2.5574373142751 2.55743731427455 2.55743731427429 2.55743731427509 2.55743731427511 2.55743731427511 2.55743731427511 2.55743731427511 2.55743731427511 2.55743731427511 2.55743731427511 2.50151548884365 2.50151548884365 2.50151548884365 2.50151548884365 2.50151548884365 2.50151548884365 2.50151548884365 2.50151548884365 2.50151548884365 2.50151548884365 2.50151548884365 2.50151548884365 2.50151548884365 2.50151548884365 2.50151548884365 2.50151548884374 2.50151548884935 2.50151548885181 2.50151548884387 2.50151548884365 2.50151548884365 2.50151548884365 2.50151548884365 2.50151548884365 2.50151548884365 2.50151548884365 2.50003135991596 2.50003135991596 2.50003135991596 2.50003135991596 2.50003135991596 2.50003135991596 2.50003135991596 2.50003135991596 2.50003135991596 2.50003135991596 2.50003135991596 2.50003135991596 2.50003135991596 2.50003135991596 2.50003135991596 2.50003135991599 2.5000313599198 2.50003135992186 2.50003135991605 2.50003135991596 2.50003135991596 2.50003135991596 2.50003135991596 2.50003135991596 2.50003135991596 2.50003135991596 2.50000052698943 2.50000052698943 2.50000052698943 2.50000052698943 2.50000052698943 2.50000052698943 2.50000052698943 2.50000052698943 2.50000052698943 2.50000052698943 2.50000052698943 2.50000052698943 2.50000052698943 2.50000052698943 2.50000052698943 2.50000052698855 2.50000052692701 2.50000052689887 2.50000052698708 2.50000052698941 2.50000052698943 2.50000052698943 2.50000052698943 2.50000052698943 2.50000052698943 2.50000052698943 2.50000000740615 2.50000000740615 2.50000000740615 2.50000000740615 2.50000000740615 2.50000000740615 2.50000000740615 2.50000000740615 2.50000000740615 2.50000000740615 2.50000000740615 2.50000000740615 2.50000000740615 2.50000000740615 2.50000000740616 2.50000000740796 2.50000000750337 2.50000000753089 2.50000000741028 2.5000000074062 2.50000000740615 2.50000000740615 2.50000000740615 2.50000000740615 2.50000000740615 2.50000000740615 2.50000000007516 2.50000000007516 2.50000000007516 2.50000000007516 2.50000000007516 2.50000000007516 2.50000000007516 2.50000000007516 2.50000000007516 2.50000000007516 2.50000000007516 2.50000000007516 2.50000000007516 2.50000000007517 2.50000000007559 2.50000000010097 2.50000000317479 2.5000000045811 2.50000000014696 2.50000000007636 2.50000000007518 2.50000000007516 2.50000000007516 2.50000000007516 2.50000000007516 2.50000000007516 2.50000000000001 2.50000000000001 2.50000000000001 2.50000000000001 2.50000000000001 2.50000000000001 2.50000000000001 2.50000000000001 2.50000000000001 2.50000000000001 2.50000000000001 2.50000000000001 2.50000000000001 2.50000000000001 2.50000000000001 2.50000000000003 2.50000000000001 2.50000000000002 2.50000000000007 2.50000000000001 2.50000000000001 2.50000000000001 2.50000000000001 2.50000000000001 2.50000000000001 2.50000000000001 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 +D/cons.5.00.000000.dat 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 +D/cons.5.00.000050.dat 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 \ No newline at end of file diff --git a/toolchain/mfc/case.py b/toolchain/mfc/case.py index 65c49f8304..7aed9212dc 100644 --- a/toolchain/mfc/case.py +++ b/toolchain/mfc/case.py @@ -295,7 +295,7 @@ def rhs_replace(match): # each element separated by new line characters. Then write those # new lines as a fully concatenated string with fortran syntax srcs.append(f"""\ - if (i == {pid}) then + if (gbl_id == {pid}) then {f"{chr(10)}".join(lines)} end if\ """) @@ -305,6 +305,7 @@ def rhs_replace(match): ! parameterize the velocity and rotation rate of a moving IB. #:def mib_analytical() +gbl_id = patch_ib(i)%gbl_patch_id {f"{chr(10)}{chr(10)}".join(srcs)} #:enddef """ diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py index 05a3fee92c..e941ffa10f 100644 --- a/toolchain/mfc/case_validator.py +++ b/toolchain/mfc/case_validator.py @@ -591,15 +591,20 @@ def check_ibm(self): ib = self.get("ib", "F") == "T" n = self.get("n", 0) num_ibs = self.get("num_ibs", 0) + num_particle_beds = self.get("num_particle_beds", 0) or 0 ib_state_wrt = self.get("ib_state_wrt", "F") == "T" self.prohibit(ib and n <= 0, "Immersed Boundaries do not work in 1D (requires n > 0)") - self.prohibit(ib and num_ibs <= 0, "num_ibs must be >= 1 when ib is enabled") - num_ib_patches_max = get_fortran_constants().get("num_ib_patches_max", 100000) + has_particle_beds = num_particle_beds > 0 and any((self.get(f"particle_bed({i})%num_particles", 0) or 0) > 0 for i in range(1, num_particle_beds + 1)) + self.prohibit( + ib and num_ibs <= 0 and not has_particle_beds, + "num_ibs must be >= 1 when ib is enabled (or specify at least one particle_bed with num_particles > 0)", + ) + num_ib_patches_max = get_fortran_constants().get("num_ib_patches_max_namelist", 50000) self.prohibit( ib and num_ibs > num_ib_patches_max, - f"num_ibs must be <= {num_ib_patches_max} (num_ib_patches_max in m_constants.fpp)", + f"num_ibs must be <= {num_ib_patches_max} (num_ib_patches_max_namelist in m_constants.fpp)", ) self.prohibit(not ib and num_ibs > 0, "num_ibs is set, but ib is not enabled") self.prohibit(ib_state_wrt and not ib, "ib_state_wrt requires ib to be enabled") diff --git a/toolchain/mfc/params/definitions.py b/toolchain/mfc/params/definitions.py index c8db5547d0..d5255b7c29 100644 --- a/toolchain/mfc/params/definitions.py +++ b/toolchain/mfc/params/definitions.py @@ -31,7 +31,8 @@ def _fc(name: str, default: int) -> int: NPR = _fc("num_probes_max", 10) # probe, acoustic, integral NB = _fc("num_bc_patches_max", 10) # patch_bc NUM_PATCHES_MAX = _fc("num_patches_max", 10) # patch_icpp (Fortran array bound) -NIB = _fc("num_ib_patches_max", 50000) # patch_ib (Fortran array bound) +NIB = _fc("num_ib_patches_max_namelist", 50000) # patch_ib namelist limit (patch_ib grows beyond this for particle beds) +NPB = _fc("num_particle_beds_max", 10) # particle_bed (Fortran array bound) # Enumeration limits for families not yet converted to IndexedFamily. # These are smaller than the Fortran array bounds to keep the registry compact. # The CONSTRAINTS dict below uses the Fortran constants for validation. @@ -229,6 +230,7 @@ def _fc(name: str, default: int) -> int: "hyperelasticity": "Enable hyperelastic model", "relativity": "Enable special relativity", "ib": "Enable immersed boundaries", + "ib_neighborhood_radius": "Neighborhood radius in ranks for IB awareness", "collision_model": "Collision model for immersed boundaries (0=none, 1=soft sphere)", "coefficient_of_restitution": "Coefficient of restitution for IB collisions", "collision_time": "Characteristic collision time for IB collisions", @@ -674,6 +676,7 @@ def get_value_label(param_name: str, value: int) -> str: "num_fluids": {"min": 1, "max": NF}, "num_patches": {"min": 0, "max": NUM_PATCHES_MAX}, "num_ibs": {"min": 0}, + "ib_neighborhood_radius": {"min": 1}, "num_source": {"min": 1}, "num_probes": {"min": 1}, "num_integrals": {"min": 1}, @@ -926,6 +929,8 @@ def _load(): # Immersed boundary _r("num_ibs", INT, {"ib"}) + _r("num_particle_beds", INT, {"ib"}) + _r("ib_neighborhood_radius", INT, {"ib"}) _r("ib", LOG, {"ib"}) _r("collision_model", INT, {"ib"}) _r("coefficient_of_restitution", REAL, {"ib"}) @@ -1179,9 +1184,8 @@ def _load(): _r(f"bub_pp%{a}", REAL, {"bubbles"}, math=sym) # patch_ib (immersed boundaries) — registered as indexed family for O(1) lookup. - # max_index is None so the parameter registry stays compact (no enumeration). - # The Fortran-side upper bound (num_ib_patches_max in m_constants.fpp) is parsed - # and enforced by the case_validator, not by max_index here. + # max_index=NIB enforces the namelist limit (num_ib_patches_max_namelist); particle beds can + # grow patch_ib beyond this at runtime, but those entries are never in the namelist. _ib_tags = {"ib"} _ib_attrs: Dict[str, tuple] = {} for a in ["geometry", "moving_ibm"]: @@ -1210,6 +1214,27 @@ def _load(): ) ) + # particle_bed — compact bed specification that expands into individual patch_ib spheres/circles at startup + _pb_tags = {"ib"} + _pb_attrs: Dict[str, tuple] = {} + for _d in ["x", "y", "z"]: + _pb_attrs[f"{_d}_centroid"] = (REAL, _pb_tags) + _pb_attrs[f"length_{_d}"] = (REAL, _pb_tags) + _pb_attrs["num_particles"] = (INT, _pb_tags) + _pb_attrs["radius"] = (REAL, _pb_tags) + _pb_attrs["mass"] = (REAL, _pb_tags) + _pb_attrs["min_spacing"] = (REAL, _pb_tags) + _pb_attrs["moving_ibm"] = (INT, _pb_tags) + _pb_attrs["seed"] = (INT, _pb_tags) + REGISTRY.register_family( + IndexedFamily( + base_name="particle_bed", + attrs=_pb_attrs, + tags=_pb_tags, + max_index=NPB, + ) + ) + # acoustic sources (4 sources) for i in range(1, NA + 1): px = f"acoustic({i})%" diff --git a/toolchain/mfc/params/descriptions.py b/toolchain/mfc/params/descriptions.py index 8ceab7f96c..952bdea0cc 100644 --- a/toolchain/mfc/params/descriptions.py +++ b/toolchain/mfc/params/descriptions.py @@ -133,6 +133,8 @@ # Immersed boundaries "ib": "Enable immersed boundary method", "num_ibs": "Number of immersed boundary patches", + "num_particle_beds": "Number of particle bed specifications to generate immersed boundary patches from", + "ib_neighborhood_radius": "Neighborhood radius in ranks for IB awareness", # Acoustic sources "acoustic_source": "Enable acoustic source terms", "num_source": "Number of acoustic sources", diff --git a/toolchain/mfc/params/namelist_parser.py b/toolchain/mfc/params/namelist_parser.py index 52385255d3..1be5b9aa78 100644 --- a/toolchain/mfc/params/namelist_parser.py +++ b/toolchain/mfc/params/namelist_parser.py @@ -208,6 +208,7 @@ "num_igr_iters", "num_igr_warm_start_iters", "num_integrals", + "num_particle_beds", "num_probes", "num_source", "nv_uvm_igr_temps_on_gpu", @@ -219,6 +220,7 @@ "p_z", "palpha_eps", "parallel_io", + "particle_bed", "patch_ib", "pi_fac", "poly_sigma", diff --git a/toolchain/mfc/run/run.py b/toolchain/mfc/run/run.py index 82e886c064..1043ab287f 100644 --- a/toolchain/mfc/run/run.py +++ b/toolchain/mfc/run/run.py @@ -52,13 +52,25 @@ def __profiler_prepend() -> typing.List[str]: return ["rocprof-compute", "profile", "-n", ARG("name").replace("-", "_").replace(".", "_")] + ARG("rcu") + ["--"] - if ARG("rsys") is not None: - if not does_command_exist("rocprof"): - raise MFCException("Failed to locate [bold red]ROCM rocprof-systems[/bold red] (rocprof-systems).") + return [] - return ["rocprof"] + ARG("rsys") - return [] +def __rsys_profiler_str() -> str: + if not does_command_exist("rocprof"): + raise MFCException("Failed to locate [bold red]ROCM rocprof-systems[/bold red] (rocprof-systems).") + + # Write a wrapper script so $SLURM_PROCID is expanded inside each srun task + # rather than by the calling shell (which would give every rank rank 0's value). + extra = shlex.join(ARG("rsys")) if ARG("rsys") else "" + wrapper_path = os.path.abspath(os.path.join(os.path.dirname(ARG("input")), "rocprof_wrapper.sh")) + wrapper_lines = [ + "#!/bin/bash", + "RANK=${SLURM_PROCID:-${FLUX_TASK_RANK:-${OMPI_COMM_WORLD_RANK:-0}}}", + f'exec rocprof -o "rocprof_rank_${{RANK}}.csv" {extra} "$@"', + ] + file_write(wrapper_path, "\n".join(wrapper_lines) + "\n") + os.chmod(wrapper_path, 0o755) + return wrapper_path def get_baked_templates() -> dict: @@ -111,7 +123,7 @@ def __generate_job_script(targets, case: input.MFCInputFile): MFC_ROOT_DIR=MFC_ROOT_DIR, SIMULATION=SIMULATION, qsystem=queues.get_system(), - profiler=shlex.join(__profiler_prepend()), + profiler=__rsys_profiler_str() if ARG("rsys") is not None else shlex.join(__profiler_prepend()), gpu_enabled=gpu_enabled, gpu_acc=gpu_acc, gpu_mp=gpu_mp, diff --git a/toolchain/templates/frontier.mako b/toolchain/templates/frontier.mako index 474baf0586..07d0a774e0 100644 --- a/toolchain/templates/frontier.mako +++ b/toolchain/templates/frontier.mako @@ -69,7 +69,10 @@ ulimit -s unlimited % if gpu_enabled: --gpus-per-task 1 --gpu-bind closest \ % endif - ${profiler} "${target.get_install_binpath(case)}") + % if target.name == 'simulation': + ${profiler} \ + % endif + "${target.get_install_binpath(case)}") % else: ${profiler} "/mnt/bb/$USER/${target.name}") % endif