diff --git a/docs/documentation/case.md b/docs/documentation/case.md index c37dbba0c0..24df0263d4 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -456,13 +456,14 @@ A particle cloud is a compact specification of a bed of identical circular (2D) | `cloud_geometry` | Integer | Shape of the cloud region. | | `shell_inner_radius` | Real | Inner radius for hemisphere-shell clouds (`cloud_geometry = 2`). | | `shell_outer_radius` | Real | Outer radius for hemisphere-shell clouds (`cloud_geometry = 2`). | +| `shell_axis` | Integer | Axis the hemisphere-shell cloud opens toward (`cloud_geometry = 2`). | | `moving_ibm` | Integer | Motion flag applied to every particle (see `patch_ib(j)%%moving_ibm`). | | `seed` | Integer | Random seed for reproducible placement (used by `packing_method = 1`). | | `packing_method` | Integer | Algorithm used to place the particles. | - `cloud_geometry` selects the cloud region: - `1` (box) uses `x[y,z]_centroid` and `length_x[y,z]` to define the region. - - `2` uses `x[y,z]_centroid`, `shell_inner_radius`, and `shell_outer_radius` to define a half-annulus in 2D and a hemisphere shell in 3D. Particle centres are sampled between `shell_inner_radius + radius` and `shell_outer_radius - radius`, and the flat plane is kept clear by one particle radius. The flat face is fixed at `y_centroid` in 2D and `z_centroid` in 3D; the filled region opens toward positive `y` in 2D and positive `z` in 3D. The full shell extent (`x[y,z]_centroid +/- shell_outer_radius` on the open side, and one particle radius of clearance on the flat-face side) must lie inside the computational domain; a hemisphere shell also requires at least two dimensions (`n > 0`). + - `2` uses `x[y,z]_centroid`, `shell_inner_radius`, `shell_outer_radius`, and `shell_axis` to define a half-annulus in 2D and a hemisphere shell in 3D. Particle centres are sampled between `shell_inner_radius + radius` and `shell_outer_radius - radius`, and the flat plane is kept clear by one particle radius. `shell_axis` (`1`=x, `2`=y, `3`=z; default `3`) selects which axis the shell opens toward from its flat face at that axis's centroid; in 2D there is no z-axis, so any value other than `1` opens toward `+y` (matching the fixed behavior before `shell_axis` existed). The open axis needs one particle radius of clearance on its flat-face side and the full `shell_outer_radius` on its open side; the other axis (2D) or two axes (3D) need the full shell extent (`centroid +/- shell_outer_radius`) inside the domain. A hemisphere shell also requires at least two dimensions (`n > 0`). - `packing_method` selects how the `num_particles` are positioned within the cloud region: - `1` (rejection sampling) draws random positions and rejects any that violate `min_spacing`, producing a disordered bed. `seed` makes the placement reproducible. - `2` (lattice) places the particles on the optimally dense lattice for the geometry — a triangular lattice in 2D and a face-centered cubic lattice in 3D. The lattice spacing is derived from the particle density (`num_particles` over the region area/volume); if that spacing is below the required `2*radius + min_spacing`, the region is too dense and the run aborts. diff --git a/src/common/m_constants.fpp b/src/common/m_constants.fpp index 40cf7d0091..09d598b19b 100644 --- a/src/common/m_constants.fpp +++ b/src/common/m_constants.fpp @@ -28,8 +28,8 @@ module m_constants integer, parameter :: num_stl_models_max = 10 !> Maximum number of immersed boundary patches (legacy, not used for patch_ib sizing) !> Fixed capacity of patch_ib (namelist patches + local particle bed subset after reduction) - integer, parameter :: num_local_ibs_max = 2000 !< Maximum number of immersed boundary patches (patch_ib) - integer, parameter :: num_ib_patches_max_namelist = 54000 + integer, parameter :: num_local_ibs_max = 8000 !< Maximum number of immersed boundary patches (patch_ib) + integer, parameter :: num_ib_patches_max_namelist = 216000 integer, parameter :: num_particle_clouds_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) diff --git a/src/common/m_derived_types.fpp b/src/common/m_derived_types.fpp index 8e66970395..12b882ab03 100644 --- a/src/common/m_derived_types.fpp +++ b/src/common/m_derived_types.fpp @@ -391,6 +391,7 @@ module m_derived_types integer :: moving_ibm !< Motion flag: 0=static, 1=moving (forces), 2=forced path integer :: seed !< Random seed for reproducible placement integer :: cloud_geometry !< Cloud region geometry: 1=box, 2=hemisphere shell + integer :: shell_axis !< Axis the hemisphere shell opens toward: 1=x, 2=y, 3=z (2D ignores 3) integer :: packing_method !< Packing algorithm: 1=rejection sampling, 2=lattice integer :: periodic !< Periodic overlap flag for box rejection packing: 0=off, 1=on end type particle_cloud_parameters @@ -520,6 +521,7 @@ module m_derived_types real(wp), dimension(3) :: ip_loc !< Physical location of the image point integer, dimension(3) :: ip_grid !< Top left grid point of IP real(wp), dimension(2, 2, 2) :: interp_coeffs !< Interpolation Coefficients of image point + logical :: interp_valid !< .false. if every image point stencil cell lies inside an IB integer :: ib_patch_id !< ID of the IB Patch the ghost point is part of real(wp) :: levelset real(wp), dimension(1:3) :: levelset_norm diff --git a/src/common/m_finite_differences.fpp b/src/common/m_finite_differences.fpp index 5639890daa..a88468b305 100644 --- a/src/common/m_finite_differences.fpp +++ b/src/common/m_finite_differences.fpp @@ -30,14 +30,29 @@ contains real(wp), dimension(-local_buff_size:q + local_buff_size), intent(in) :: s_cc integer :: i !< Generic loop iterator + ! Coefficients always extend at least fd_number_in beyond the interior on each side, so a stencil centered on a + ! ghost-adjacent cell (e.g. an immersed boundary near a domain boundary) has a real coefficient to read instead of + ! reading past the caller's allocation. offset_s, when given, widens this further (never narrows it) for callers + ! that need more than fd_number_in of margin. + if (present(offset_s)) then - lB = -offset_s%beg - lE = q + offset_s%end + lB = -max(fd_number_in, offset_s%beg) + lE = q + max(fd_number_in, offset_s%end) else - lB = 0 - lE = q + lB = -fd_number_in + lE = q + fd_number_in end if + ! The stencil below reaches fd_number_in cells either side of i, and s_cc only exists over the + ! caller's buffer, so the loop cannot start closer than fd_number_in to either end of it. Without + ! this, a caller whose buffer is narrower than 2*fd_number_in reads off the front of s_cc: a + ! 4th-order hypoelastic case with buff_size = 2 starts at i = -2 and immediately asks for + ! s_cc(-4). Immersed-boundary cases are unaffected -- s_mfc_buff_size gives them buff_size >= 10, + ! so the clamp never binds and they still get coefficients out to -fd_number_in, which is the + ! whole point of widening the range above. + lB = max(lB, -local_buff_size + fd_number_in) + lE = min(lE, q + local_buff_size - fd_number_in) + ! Computing the 1st order finite-difference coefficients if (fd_order_in == 1) then do i = lB, lE diff --git a/src/post_process/m_data_output.fpp b/src/post_process/m_data_output.fpp index 795ae7e736..8151325a19 100644 --- a/src/post_process/m_data_output.fpp +++ b/src/post_process/m_data_output.fpp @@ -454,31 +454,54 @@ contains ! Finally, the local quadrilateral mesh, either 2D or 3D, along with its offsets that indicate the presence and size of ! ghost zone layer(s), are put in the formatted database slave file. + ! Silo carries the mesh coordinates in their own datatype, separate from the flow variables, so the cell boundaries + ! are copied down to single precision here when that is what was asked for. Without this the mesh is always written as + ! DB_DOUBLE, which keeps downstream readers on a double-precision path regardless of `precision`. + if (precision == precision_single) then + x_cb_s = real(x_cb, sp) + if (n > 0) then + y_cb_s = real(y_cb, sp) + if (p > 0) z_cb_s = real(z_cb, sp) + end if + end if + if (p > 0) then err = DBMKOPTLIST(2, out%optlist) err = DBADDIAOPT(out%optlist, DBOPT_LO_OFFSET, size(out%lo_offset), out%lo_offset) err = DBADDIAOPT(out%optlist, DBOPT_HI_OFFSET, size(out%hi_offset), out%hi_offset) - if (grid_geometry == 3) then - err = DBPUTQM(out%dbfile, 'rectilinear_grid', 16, 'x', 1, 'y', 1, 'z', 1, y_cb, z_cb, x_cb, out%dims, 3, & - & DB_DOUBLE, DB_COLLINEAR, out%optlist, ierr) - else - err = DBPUTQM(out%dbfile, 'rectilinear_grid', 16, 'x', 1, 'y', 1, 'z', 1, x_cb, y_cb, z_cb, out%dims, 3, & - & DB_DOUBLE, DB_COLLINEAR, out%optlist, ierr) - end if + #:for PRECISION, SFX, DBT in [(1,'_s','DB_FLOAT'),(2,'',"DB_DOUBLE")] + if (precision == ${PRECISION}$) then + if (grid_geometry == 3) then + err = DBPUTQM(out%dbfile, 'rectilinear_grid', 16, 'x', 1, 'y', 1, 'z', 1, y_cb${SFX}$, z_cb${SFX}$, & + & x_cb${SFX}$, out%dims, 3, ${DBT}$, DB_COLLINEAR, out%optlist, ierr) + else + err = DBPUTQM(out%dbfile, 'rectilinear_grid', 16, 'x', 1, 'y', 1, 'z', 1, x_cb${SFX}$, y_cb${SFX}$, & + & z_cb${SFX}$, out%dims, 3, ${DBT}$, DB_COLLINEAR, out%optlist, ierr) + end if + end if + #:endfor err = DBFREEOPTLIST(out%optlist) else if (n > 0) then err = DBMKOPTLIST(2, out%optlist) err = DBADDIAOPT(out%optlist, DBOPT_LO_OFFSET, size(out%lo_offset), out%lo_offset) err = DBADDIAOPT(out%optlist, DBOPT_HI_OFFSET, size(out%hi_offset), out%hi_offset) - err = DBPUTQM(out%dbfile, 'rectilinear_grid', 16, 'x', 1, 'y', 1, 'z', 1, x_cb, y_cb, DB_F77NULL, out%dims, 2, & - & DB_DOUBLE, DB_COLLINEAR, out%optlist, ierr) + #:for PRECISION, SFX, DBT in [(1,'_s','DB_FLOAT'),(2,'',"DB_DOUBLE")] + if (precision == ${PRECISION}$) then + err = DBPUTQM(out%dbfile, 'rectilinear_grid', 16, 'x', 1, 'y', 1, 'z', 1, x_cb${SFX}$, y_cb${SFX}$, & + & DB_F77NULL, out%dims, 2, ${DBT}$, DB_COLLINEAR, out%optlist, ierr) + end if + #:endfor err = DBFREEOPTLIST(out%optlist) else err = DBMKOPTLIST(2, out%optlist) err = DBADDIAOPT(out%optlist, DBOPT_LO_OFFSET, size(out%lo_offset), out%lo_offset) err = DBADDIAOPT(out%optlist, DBOPT_HI_OFFSET, size(out%hi_offset), out%hi_offset) - err = DBPUTQM(out%dbfile, 'rectilinear_grid', 16, 'x', 1, 'y', 1, 'z', 1, x_cb, DB_F77NULL, DB_F77NULL, out%dims, & - & 1, DB_DOUBLE, DB_COLLINEAR, out%optlist, ierr) + #:for PRECISION, SFX, DBT in [(1,'_s','DB_FLOAT'),(2,'',"DB_DOUBLE")] + if (precision == ${PRECISION}$) then + err = DBPUTQM(out%dbfile, 'rectilinear_grid', 16, 'x', 1, 'y', 1, 'z', 1, x_cb${SFX}$, DB_F77NULL, & + & DB_F77NULL, out%dims, 1, ${DBT}$, DB_COLLINEAR, out%optlist, ierr) + end if + #:endfor err = DBFREEOPTLIST(out%optlist) end if else if (format == format_binary) then @@ -1349,6 +1372,10 @@ contains real(wp), dimension(:), allocatable :: omega_x, omega_y, omega_z real(wp), dimension(:), allocatable :: angle_x, angle_y, angle_z real(wp), dimension(:), allocatable :: ib_diameter + real(sp), dimension(:), allocatable :: px_s, py_s, pz_s + logical, dimension(:), allocatable :: keep + integer :: nKept + real(wp) :: r_ib if (proc_rank == 0) then nBodies = num_ibs @@ -1419,12 +1446,51 @@ contains ib_diameter(i) = ib_data(i, 20)*2.0_wp end do + ! When only part of the domain is written, the bodies outside that window are dropped so the point mesh matches + ! the cropped grid. A body is kept when its bounding sphere overlaps the window rather than when its centroid is + ! inside it, so one straddling the boundary still appears instead of vanishing at the edge. + if (output_partial_domain) then + allocate (keep(nBodies)) + + do i = 1, nBodies + r_ib = 0.5_wp*ib_diameter(i) + keep(i) = (px(i) + r_ib >= x_output%beg) .and. (px(i) - r_ib <= x_output%end) + if (n > 0) keep(i) = keep(i) .and. (py(i) + r_ib >= y_output%beg) .and. (py(i) - r_ib <= y_output%end) + if (p > 0) keep(i) = keep(i) .and. (pz(i) + r_ib >= z_output%beg) .and. (pz(i) - r_ib <= z_output%end) + end do + + nKept = count(keep) + + if (nKept < nBodies) then + #:for A in ['px','py','pz','ib_diameter','force_x','force_y','force_z','torque_x','torque_y','torque_z'] + ${A}$(1:nKept) = pack(${A}$(1:nBodies), keep) + #:endfor + #:for A in ['vel_x','vel_y','vel_z','omega_x','omega_y','omega_z','angle_x','angle_y','angle_z'] + ${A}$(1:nKept) = pack(${A}$(1:nBodies), keep) + #:endfor + end if + + nBodies = nKept + deallocate (keep) + end if + write (meshnames(1), '(A,I0,A)') '../p0/', t_step, '.silo:ib_bodies' meshtypes(1) = DB_POINTMESH err = DBSET2DSTRLEN(len(meshnames(1))) err = DBPUTMMESH(out%dbroot, 'ib_bodies', 16, 1, meshnames, len_trim(meshnames), meshtypes, DB_F77NULL, ierr) - err = DBPUTPM(out%dbfile, 'ib_bodies', 9, 3, px, py, pz, nBodies, DB_DOUBLE, DB_F77NULL, ierr) + ! Silo carries the point-mesh coordinates in their own datatype, so they need the same single-precision + ! treatment as the rectilinear mesh. + if (precision == precision_single) then + allocate (px_s(nBodies), py_s(nBodies), pz_s(nBodies)) + px_s(1:nBodies) = real(px(1:nBodies), sp) + py_s(1:nBodies) = real(py(1:nBodies), sp) + pz_s(1:nBodies) = real(pz(1:nBodies), sp) + err = DBPUTPM(out%dbfile, 'ib_bodies', 9, 3, px_s, py_s, pz_s, nBodies, DB_FLOAT, DB_F77NULL, ierr) + deallocate (px_s, py_s, pz_s) + else + err = DBPUTPM(out%dbfile, 'ib_bodies', 9, 3, px, py, pz, nBodies, DB_DOUBLE, DB_F77NULL, ierr) + end if 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) @@ -1456,12 +1522,13 @@ 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) :: var_name_entry - integer :: var_type_entry, ierr + 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 + real(sp), dimension(:), allocatable :: data_s write (var_name_entry, '(A,I0,A)') '../p0/', t_step, '.silo:' // trim(varname) var_type_entry = DB_POINTVAR @@ -1469,7 +1536,15 @@ contains err = DBPUTMVAR(out%dbroot, trim(varname), len_trim(varname), 1, var_name_entry, len_trim(var_name_entry), & & var_type_entry, DB_F77NULL, ierr) - err = DBPUTPV1(out%dbfile, trim(varname), len_trim(varname), 'ib_bodies', 9, data, nBodies, DB_DOUBLE, DB_F77NULL, ierr) + if (precision == precision_single) then + allocate (data_s(nBodies)) + data_s(1:nBodies) = real(data(1:nBodies), sp) + err = DBPUTPV1(out%dbfile, trim(varname), len_trim(varname), 'ib_bodies', 9, data_s, nBodies, DB_FLOAT, DB_F77NULL, & + & ierr) + deallocate (data_s) + else + err = DBPUTPV1(out%dbfile, trim(varname), len_trim(varname), 'ib_bodies', 9, data, nBodies, DB_DOUBLE, DB_F77NULL, ierr) + end if end subroutine s_write_ib_variable diff --git a/src/post_process/m_derived_variables.fpp b/src/post_process/m_derived_variables.fpp index 105512143e..6d5e044c96 100644 --- a/src/post_process/m_derived_variables.fpp +++ b/src/post_process/m_derived_variables.fpp @@ -33,18 +33,20 @@ contains allocate (fd%gm_rho_sf(-offset_x%beg:m + offset_x%end,-offset_y%beg:n + offset_y%end,-offset_z%beg:p + offset_z%end)) end if - ! Allocate FD coefficients (up to 4th order; higher orders need extension) + ! Allocate FD coefficients (up to 4th order; higher orders need extension). s_compute_finite_difference_coefficients + ! always extends at least fd_number beyond the interior on each side, widened further by offset_x/y/z when those + ! are larger (multi-block Silo ghost zones); the allocation must cover whichever bound ends up wider. if (omega_wrt(2) .or. omega_wrt(3) .or. qm_wrt .or. schlieren_wrt .or. liutex_wrt) then - allocate (fd%fd_coeff_x(-fd_number:fd_number,-offset_x%beg:m + offset_x%end)) + allocate (fd%fd_coeff_x(-fd_number:fd_number,-max(fd_number, offset_x%beg):m + max(fd_number, offset_x%end))) end if if (omega_wrt(1) .or. omega_wrt(3) .or. qm_wrt .or. liutex_wrt .or. (n > 0 .and. schlieren_wrt)) then - allocate (fd%fd_coeff_y(-fd_number:fd_number,-offset_y%beg:n + offset_y%end)) + allocate (fd%fd_coeff_y(-fd_number:fd_number,-max(fd_number, offset_y%beg):n + max(fd_number, offset_y%end))) end if if (omega_wrt(1) .or. omega_wrt(2) .or. qm_wrt .or. liutex_wrt .or. (p > 0 .and. schlieren_wrt)) then - allocate (fd%fd_coeff_z(-fd_number:fd_number,-offset_z%beg:p + offset_z%end)) + allocate (fd%fd_coeff_z(-fd_number:fd_number,-max(fd_number, offset_z%beg):p + max(fd_number, offset_z%end))) end if end subroutine s_initialize_derived_variables_module diff --git a/src/post_process/m_global_parameters.fpp b/src/post_process/m_global_parameters.fpp index 2a8c9e9929..bd90c28ed8 100644 --- a/src/post_process/m_global_parameters.fpp +++ b/src/post_process/m_global_parameters.fpp @@ -51,6 +51,11 @@ module m_global_parameters !> @name Cell-boundary locations in the x-, y- and z-coordinate directions !> @{ real(wp), allocatable, dimension(:) :: x_cb, x_root_cb, y_cb, z_cb + ! Single-precision copies, handed to Silo when precision == precision_single. + ! Silo stores the mesh coordinates with their own datatype, independent of the + ! flow variables, so the mesh needs its own single-precision arrays to follow + ! the requested precision. + real(sp), allocatable, dimension(:) :: x_cb_s, y_cb_s, z_cb_s !> @} !> @name Cell-center locations in the x-, y- and z-coordinate directions @@ -521,16 +526,28 @@ contains allocate (x_cc(-buff_size:m + buff_size)) allocate (dx(-buff_size:m + buff_size)) + if (precision == precision_single) then + allocate (x_cb_s(-1 - offset_x%beg:m + offset_x%end)) + end if + ! Allocating grid variables in the y- and z-coordinate directions if (n > 0) then allocate (y_cb(-1 - offset_y%beg:n + offset_y%end)) allocate (y_cc(-buff_size:n + buff_size)) allocate (dy(-buff_size:n + buff_size)) + if (precision == precision_single) then + allocate (y_cb_s(-1 - offset_y%beg:n + offset_y%end)) + end if + if (p > 0) then allocate (z_cb(-1 - offset_z%beg:p + offset_z%end)) allocate (z_cc(-buff_size:p + buff_size)) allocate (dz(-buff_size:p + buff_size)) + + if (precision == precision_single) then + allocate (z_cb_s(-1 - offset_z%beg:p + offset_z%end)) + end if end if ! Allocating the grid variables, only used for the 1D simulations, and containing the defragmented computational domain @@ -575,12 +592,15 @@ contains ! Deallocating the grid variables for the x-coordinate direction deallocate (x_cc, x_cb, dx) + if (allocated(x_cb_s)) deallocate (x_cb_s) ! Deallocating grid variables for the y- and z-coordinate directions if (n > 0) then deallocate (y_cc, y_cb, dy) + if (allocated(y_cb_s)) deallocate (y_cb_s) if (p > 0) then deallocate (z_cc, z_cb, dz) + if (allocated(z_cb_s)) deallocate (z_cb_s) end if else ! Deallocating the grid variables, only used for the 1D simulations, and containing the defragmented computational diff --git a/src/pre_process/m_particle_cloud.fpp b/src/pre_process/m_particle_cloud.fpp index fba6452ca9..6119079fb5 100644 --- a/src/pre_process/m_particle_cloud.fpp +++ b/src/pre_process/m_particle_cloud.fpp @@ -22,7 +22,7 @@ module m_particle_cloud private - public :: s_generate_particle_clouds + public :: s_generate_particle_clouds, s_add_cloud_particle contains @@ -182,8 +182,9 @@ contains !> Draws one rejection-sampling candidate centre (rx, ry, rz) for cloud_idx, advancing seed in place. For box geometry the !! candidate is uniform in the box and never rejected. For a hemisphere shell the candidate is uniform in the shell volume - 2D !! uses theta uniform on [0, pi] with the sqrt radial CDF; 3D uses uniform phi, uniform cos(polar) on [0, 1], and the cube-root - !! radial CDF - and reject is set when it lands within one particle radius of the flat face (the plane at y_centroid in 2D, - !! z_centroid in 3D), a hard geometric cut applied after sampling that preserves uniformity over the remaining region. + !! radial CDF - and reject is set when it lands within one particle radius of the flat face (the plane through the centroid + !! perpendicular to shell_axis: 1=x, 2=y, 3=z; 2D has no z-axis so any value other than 1 falls back to y, matching the + !! pre-shell_axis default), a hard geometric cut applied after sampling that preserves uniformity over the remaining region. subroutine s_sample_cloud_candidate(cloud_idx, seed, rx, ry, rz, reject) integer, intent(in) :: cloud_idx @@ -219,20 +220,39 @@ contains theta = pi*f_xorshift(seed) u = f_xorshift(seed) r_shell = sqrt((r_outer**2 - r_inner**2)*u + r_inner**2) - rx = particle_cloud(cloud_idx)%x_centroid + r_shell*cos(theta) - ry = particle_cloud(cloud_idx)%y_centroid + r_shell*sin(theta) rz = particle_cloud(cloud_idx)%z_centroid - if (ry < particle_cloud(cloud_idx)%y_centroid + particle_cloud(cloud_idx)%radius) reject = .true. + if (particle_cloud(cloud_idx)%shell_axis == 1) then + rx = particle_cloud(cloud_idx)%x_centroid + r_shell*sin(theta) + ry = particle_cloud(cloud_idx)%y_centroid + r_shell*cos(theta) + if (rx < particle_cloud(cloud_idx)%x_centroid + particle_cloud(cloud_idx)%radius) reject = .true. + else + rx = particle_cloud(cloud_idx)%x_centroid + r_shell*cos(theta) + ry = particle_cloud(cloud_idx)%y_centroid + r_shell*sin(theta) + if (ry < particle_cloud(cloud_idx)%y_centroid + particle_cloud(cloud_idx)%radius) reject = .true. + end if else phi = 2._wp*pi*f_xorshift(seed) zdir = f_xorshift(seed) rho = sqrt(max(0._wp, 1._wp - zdir**2)) u = f_xorshift(seed) r_shell = ((r_outer**3 - r_inner**3)*u + r_inner**3)**(1._wp/3._wp) - rx = particle_cloud(cloud_idx)%x_centroid + r_shell*rho*cos(phi) - ry = particle_cloud(cloud_idx)%y_centroid + r_shell*rho*sin(phi) - rz = particle_cloud(cloud_idx)%z_centroid + r_shell*zdir - if (rz < particle_cloud(cloud_idx)%z_centroid + particle_cloud(cloud_idx)%radius) reject = .true. + select case (particle_cloud(cloud_idx)%shell_axis) + case (1) ! opens toward +x + rx = particle_cloud(cloud_idx)%x_centroid + r_shell*zdir + ry = particle_cloud(cloud_idx)%y_centroid + r_shell*rho*cos(phi) + rz = particle_cloud(cloud_idx)%z_centroid + r_shell*rho*sin(phi) + if (rx < particle_cloud(cloud_idx)%x_centroid + particle_cloud(cloud_idx)%radius) reject = .true. + case (2) ! opens toward +y + ry = particle_cloud(cloud_idx)%y_centroid + r_shell*zdir + rx = particle_cloud(cloud_idx)%x_centroid + r_shell*rho*cos(phi) + rz = particle_cloud(cloud_idx)%z_centroid + r_shell*rho*sin(phi) + if (ry < particle_cloud(cloud_idx)%y_centroid + particle_cloud(cloud_idx)%radius) reject = .true. + case default ! 3: opens toward +z + rz = particle_cloud(cloud_idx)%z_centroid + r_shell*zdir + rx = particle_cloud(cloud_idx)%x_centroid + r_shell*rho*cos(phi) + ry = particle_cloud(cloud_idx)%y_centroid + r_shell*rho*sin(phi) + if (rz < particle_cloud(cloud_idx)%z_centroid + particle_cloud(cloud_idx)%radius) reject = .true. + end select end if case default call s_mpi_abort("Particle cloud geometry is not a known cloud geometry of MFC. Exiting.") diff --git a/src/simulation/m_bubbles_EL.fpp b/src/simulation/m_bubbles_EL.fpp index 15a709c179..35b1b3f812 100644 --- a/src/simulation/m_bubbles_EL.fpp +++ b/src/simulation/m_bubbles_EL.fpp @@ -186,18 +186,19 @@ contains if (lag_params%vel_model > 0 .and. lag_params%pressure_force) then @:ALLOCATE(grad_p_x(0:m, 0:n, 0:p)) - @:ALLOCATE(fd_coeff_x_pgrad(-fd_number:fd_number, 0:m)) + ! s_compute_finite_difference_coefficients always extends fd_number beyond the interior on each side + @:ALLOCATE(fd_coeff_x_pgrad(-fd_number:fd_number,-fd_number:m + fd_number)) call s_compute_finite_difference_coefficients(m, x_cc, fd_coeff_x_pgrad, buff_size, fd_number, fd_order) $:GPU_UPDATE(device='[fd_coeff_x_pgrad]') if (n > 0) then @:ALLOCATE(grad_p_y(0:m, 0:n, 0:p)) - @:ALLOCATE(fd_coeff_y_pgrad(-fd_number:fd_number, 0:n)) + @:ALLOCATE(fd_coeff_y_pgrad(-fd_number:fd_number,-fd_number:n + fd_number)) call s_compute_finite_difference_coefficients(n, y_cc, fd_coeff_y_pgrad, buff_size, fd_number, fd_order) $:GPU_UPDATE(device='[fd_coeff_y_pgrad]') end if if (p > 0) then @:ALLOCATE(grad_p_z(0:m, 0:n, 0:p)) - @:ALLOCATE(fd_coeff_z_pgrad(-fd_number:fd_number, 0:p)) + @:ALLOCATE(fd_coeff_z_pgrad(-fd_number:fd_number,-fd_number:p + fd_number)) call s_compute_finite_difference_coefficients(p, z_cc, fd_coeff_z_pgrad, buff_size, fd_number, fd_order) $:GPU_UPDATE(device='[fd_coeff_z_pgrad]') end if diff --git a/src/simulation/m_collisions.fpp b/src/simulation/m_collisions.fpp index c8148fcd15..1a7a24f589 100644 --- a/src/simulation/m_collisions.fpp +++ b/src/simulation/m_collisions.fpp @@ -91,7 +91,7 @@ contains real(wp), dimension(num_ibs, 3), intent(inout) :: forces, torques integer :: i, encoded_pid1, encoded_pid2, xp1, xp2, yp1, yp2, zp1, zp2, pid1, pid2, l ! iterators and patch IDs real(wp) :: overlap_distance - real(wp), dimension(3) :: normal_vector, centroid_1, centroid_2 + real(wp), dimension(3) :: normal_vector, centroid_1, centroid_2, contact_point real(wp), dimension(3) :: normal_velocity, tangential_vector, normal_force, tangential_force, torque, radial_vector, & & rotation_velocity, vel1, vel2 real(wp) :: k, eta, effective_mass ! the spring stiffness and damping coefficient and mass of a specific interaction @@ -102,7 +102,7 @@ contains $:GPU_PARALLEL_LOOP(private='[i, l, encoded_pid1, encoded_pid2, xp1, xp2, yp1, yp2, zp1, zp2, pid1, pid2, centroid_1, & & centroid_2, normal_vector, overlap_distance, effective_mass, k, eta, normal_velocity, & & tangential_vector, normal_force, tangential_force, torque, radial_vector, rotation_velocity, vel1, & - & vel2]', copy='[forces, torques]') + & vel2, contact_point]', copy='[forces, torques]') do i = 1, num_considered_collisions encoded_pid1 = collision_lookup(i, 3) encoded_pid2 = collision_lookup(i, 4) @@ -128,7 +128,9 @@ 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_location(centroid_1, glb_bounds)) then + ! pid1 is a rank-local index, so owning the pair by its centroid drops or doubles pairs split across ranks + contact_point = centroid_1 + normal_vector*(patch_ib(pid1)%radius - 0.5_wp*overlap_distance) + if (f_local_rank_owns_location(contact_point, glb_bounds)) 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 @@ -274,14 +276,16 @@ contains do kk = k - z_bound, k + z_bound neighbor_patch_id = ib_markers%sf(ii, jj, kk) - ! If any neighbors are of a different/higher marker value, we consider it for possible collision - if (gp_patch_id < neighbor_patch_id) then + ! Any neighbor of a different patch is a candidate pair. Both patches record it: the rank that owns the + ! contact point may hold interior ghost points of only one of them, so one-sided detection can leave + ! that rank blind to a contact that sits within a cell of its boundary. The host pass below sorts the + ! pair and drops duplicates. + if (neighbor_patch_id /= 0 .and. neighbor_patch_id /= gp_patch_id) then $:GPU_ATOMIC(atomic='capture') num_raw = num_raw + 1 local_num_raw = num_raw $:END_GPU_ATOMIC_CAPTURE() - ! Store with smaller ID first for consistent ordering raw_pairs(local_num_raw, 1) = gp_patch_id raw_pairs(local_num_raw, 2) = neighbor_patch_id exit neighbor_search diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index dbbb4445c6..a1f3425804 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -185,6 +185,7 @@ contains real(wp) :: icfl, vcfl, ccfl, tcfl, Rc real(wp) :: mu_frac, mu_frac_max_loc, mu_frac_max_glb !< Compression as a fraction of the EOS limit integer :: fl !< Fluid loop iterator + logical :: include_cell !< Cell is fluid, not ghost/inside an IB real(wp), dimension(4) :: stab_max_loc, stab_max_glb !< Max-reduced criteria (ICFL, VCFL, CCFL, TCFL), packed real(wp), dimension(1) :: stab_min_loc, stab_min_glb !< Min-reduced criteria (Rc), packed @@ -196,47 +197,54 @@ contains mu_frac_max_loc = 0._wp ! Computing Stability Criteria at Current Time-step $:GPU_PARALLEL_LOOP(collapse=3, private='[j, k, l, vel, alpha, alpha_rho, Re, rho, vel_sum, pres, gamma, pi_inf, c, qv, & - & icfl, vcfl, Rc, ccfl, tcfl, fl, mu_frac]', reduction='[[icfl_max_loc, vcfl_max_loc, ccfl_max_loc, & - & tcfl_max_loc, mu_frac_max_loc], [Rc_min_loc]]', reductionOp='[max, min]') + & icfl, vcfl, Rc, ccfl, tcfl, fl, mu_frac, include_cell]', reduction='[[icfl_max_loc, vcfl_max_loc, & + & ccfl_max_loc, tcfl_max_loc, mu_frac_max_loc], [Rc_min_loc]]', reductionOp='[max, min]') do l = 0, p do k = 0, n do j = 0, m - call s_compute_cell_state(q_prim_vf, pres, rho, gamma, pi_inf, Re, alpha, alpha_rho, vel, vel_sum, qv, j, k, l) + ! exclude cells inside of immersed boundaries + include_cell = .true. + if (ib) include_cell = (ib_markers%sf(j, k, l) == 0) + if (include_cell) then + call s_compute_cell_state(q_prim_vf, pres, rho, gamma, pi_inf, Re, alpha, alpha_rho, vel, vel_sum, qv, j, & + & k, l) - call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, alpha, c, alpha_rho) + call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, alpha, c, alpha_rho) - ! How close each Mie-Gruneisen phase is to the compression its Hugoniot fit can represent. - ! Past 1 there is no shock state to find and the reference curve is fiction, so it is reduced - ! out of the kernel and turned into an abort on the host -- s_mpi_abort cannot be called here. - if (any_state_dependent_eos) then - $:GPU_LOOP(parallelism='[seq]') - do fl = 1, num_fluids - if (eoss(fl) == eos_mie_gruneisen) then - mu_frac = (alpha_rho(fl)/max(alpha(fl), sgm_eps)/eos_coeffs(fl)%rho0 - 1._wp)/eos_coeffs(fl)%mu_max - mu_frac_max_loc = max(mu_frac_max_loc, mu_frac) - end if - end do - end if + ! How close each Mie-Gruneisen phase is to the compression its Hugoniot fit can represent. + ! Past 1 there is no shock state to find and the reference curve is fiction, so it is reduced + ! out of the kernel and turned into an abort on the host -- s_mpi_abort cannot be called here. + if (any_state_dependent_eos) then + $:GPU_LOOP(parallelism='[seq]') + do fl = 1, num_fluids + if (eoss(fl) == eos_mie_gruneisen) then + mu_frac = (alpha_rho(fl)/max(alpha(fl), & + & sgm_eps)/eos_coeffs(fl)%rho0 - 1._wp)/eos_coeffs(fl)%mu_max + mu_frac_max_loc = max(mu_frac_max_loc, mu_frac) + end if + end do + end if - if (any_non_newtonian) then - Re(1) = 0._wp - do fl = 1, num_fluids - if (is_non_newtonian(fl)) then - Re(1) = Re(1) + alpha(fl)*hb_mu_max(fl) - else - Re(1) = Re(1) + alpha(fl)*fluid_inv_re(fl) - end if - end do - Re(1) = 1._wp/max(Re(1), sgm_eps) - end if + if (any_non_newtonian) then + Re(1) = 0._wp + do fl = 1, num_fluids + if (is_non_newtonian(fl)) then + Re(1) = Re(1) + alpha(fl)*hb_mu_max(fl) + else + Re(1) = Re(1) + alpha(fl)*fluid_inv_re(fl) + end if + end do + Re(1) = 1._wp/max(Re(1), sgm_eps) + end if - call s_compute_stability_from_dt(vel, c, rho, Re, alpha, alpha_rho, j, k, l, icfl, vcfl, Rc, ccfl, tcfl) + call s_compute_stability_from_dt(vel, c, rho, Re, alpha, alpha_rho, j, k, l, icfl, vcfl, Rc, ccfl, tcfl) - icfl_max_loc = max(icfl_max_loc, icfl) - vcfl_max_loc = max(vcfl_max_loc, merge(vcfl, 0.0_wp, viscous)) - ccfl_max_loc = max(ccfl_max_loc, merge(ccfl, 0.0_wp, surface_tension)) - tcfl_max_loc = max(tcfl_max_loc, merge(tcfl, 0.0_wp, heat_conduction)) - Rc_min_loc = min(Rc_min_loc, merge(Rc, huge(1.0_wp), viscous)) + icfl_max_loc = max(icfl_max_loc, icfl) + vcfl_max_loc = max(vcfl_max_loc, merge(vcfl, 0.0_wp, viscous)) + ccfl_max_loc = max(ccfl_max_loc, merge(ccfl, 0.0_wp, surface_tension)) + tcfl_max_loc = max(tcfl_max_loc, merge(tcfl, 0.0_wp, heat_conduction)) + Rc_min_loc = min(Rc_min_loc, merge(Rc, huge(1.0_wp), viscous)) + end if end do end do end do @@ -280,6 +288,13 @@ contains if (Rc_min_glb < Rc_min) Rc_min = Rc_min_glb end if + ! Any rank whose own local extremum violates the limit is, by construction of the + ! max-reduction above, a rank that actually contains the offending cell(s). + if ((.not. f_approx_equal(icfl_max_loc, icfl_max_loc)) .or. icfl_max_loc > 1._wp) then + call s_report_icfl_violation(q_prim_vf) + end if + call s_mpi_barrier() ! ensure diagnostic output above is flushed before any rank aborts below + if (proc_rank == 0) then write (3, '(13X,I9,13X,F10.6,13X,F10.6,13X,F10.6)', advance="no") t_step, dt, mytime, icfl_max_glb @@ -333,6 +348,135 @@ contains end subroutine s_write_run_time_information + !> Locate the grid cell responsible for an ICFL violation on this rank and report its state plus the nearest immersed-boundary + !! particles, to aid debugging stability failures in particle-laden high-Mach cases. + impure subroutine s_report_icfl_violation(q_prim_vf) + + type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf + real(wp), dimension(num_fluids) :: alpha, alpha_rho + real(wp), dimension(num_vels) :: vel, vel_hit + real(wp), dimension(2) :: Re + real(wp) :: rho, vel_sum, pres, gamma, pi_inf, qv, c + real(wp) :: rho_hit, pres_hit, c_hit + real(wp) :: icfl, vcfl, Rc, ccfl, tcfl, icfl_hit + integer :: i, j, k, l, fl, j_hit, k_hit, l_hit + real(wp) :: x_hit, y_hit, z_hit, dist + logical :: nan_hit + integer :: near1_id, near2_id + real(wp) :: near1_dist, near2_dist + + do i = 1, sys_size + $:GPU_UPDATE(host='[q_prim_vf(i)%sf(:, :, :)]') + end do + if (ib) then + $:GPU_UPDATE(host='[ib_markers%sf]') + end if + + icfl_hit = -huge(1._wp) + nan_hit = .false. + j_hit = 0; k_hit = 0; l_hit = 0 + + scan: do l = 0, p + do k = 0, n + do j = 0, m + if (ib) then + if (ib_markers%sf(j, k, l) /= 0) cycle + end if + + call s_compute_cell_state(q_prim_vf, pres, rho, gamma, pi_inf, Re, alpha, alpha_rho, vel, vel_sum, qv, j, k, l) + call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, alpha, c, alpha_rho) + + if (any_non_newtonian) then + Re(1) = 0._wp + do fl = 1, num_fluids + if (is_non_newtonian(fl)) then + Re(1) = Re(1) + alpha(fl)*hb_mu_max(fl) + else + Re(1) = Re(1) + alpha(fl)*fluid_inv_re(fl) + end if + end do + Re(1) = 1._wp/max(Re(1), sgm_eps) + end if + + call s_compute_stability_from_dt(vel, c, rho, Re, alpha, alpha_rho, j, k, l, icfl, vcfl, Rc, ccfl, tcfl) + + if (.not. f_approx_equal(icfl, icfl)) then + nan_hit = .true. + j_hit = j; k_hit = k; l_hit = l + rho_hit = rho; pres_hit = pres; c_hit = c; vel_hit = vel + exit scan + else if (icfl > icfl_hit) then + icfl_hit = icfl + j_hit = j; k_hit = k; l_hit = l + rho_hit = rho; pres_hit = pres; c_hit = c; vel_hit = vel + end if + end do + end do + end do scan + + x_hit = x_cc(j_hit) + y_hit = 0._wp; if (n > 0) y_hit = y_cc(k_hit) + z_hit = 0._wp; if (p > 0) z_hit = z_cc(l_hit) + + print '(A,I0,A,I0,A,I0,A,I0,A)', 'ICFL violation on rank ', proc_rank, ': cell (j,k,l) = (', j_hit, ',', k_hit, ',', & + & l_hit, ')' + if (nan_hit) then + print '(A)', ' icfl = NaN' + else + print '(A,ES16.6)', ' icfl = ', icfl_hit + end if + print '(A,3(ES16.6,1X))', ' position = ', x_hit, y_hit, z_hit + print '(A,ES16.6,A,ES16.6,A,ES16.6)', ' rho, pres, c = ', rho_hit, ', ', pres_hit, ', ', c_hit + print '(A,3(ES16.6,1X))', ' velocity = ', vel_hit + if (ib) print '(A,I0)', ' ib_markers = ', ib_markers%sf(j_hit, k_hit, l_hit) + + if (ib .and. num_ibs > 0) then + near1_id = 0; near1_dist = huge(1._wp) + near2_id = 0; near2_dist = huge(1._wp) + do i = 1, num_ibs + dist = sqrt((x_hit - patch_ib(i)%x_centroid)**2 + (y_hit - patch_ib(i)%y_centroid)**2 + (z_hit & + & - patch_ib(i)%z_centroid)**2) + if (dist < near1_dist) then + near2_dist = near1_dist; near2_id = near1_id + near1_dist = dist; near1_id = i + else if (dist < near2_dist) then + near2_dist = dist; near2_id = i + end if + end do + if (near1_id > 0) then + print '(A,I0,A,ES16.6,A,ES16.6,A,3(ES16.6,1X))', ' nearest particle id=', near1_id, ' dist=', near1_dist, & + & ' gap=', near1_dist - patch_ib(near1_id)%radius, ' vel=', patch_ib(near1_id)%vel + print '(A,3(ES16.6,1X))', ' centroid = ', patch_ib(near1_id)%x_centroid, patch_ib(near1_id)%y_centroid, & + & patch_ib(near1_id)%z_centroid + print '(A,3(ES16.6,1X))', ' angular_vel = ', patch_ib(near1_id)%angular_vel + print '(A,3(ES16.6,1X))', ' force = ', patch_ib(near1_id)%force + print '(A,3(ES16.6,1X))', ' torque = ', patch_ib(near1_id)%torque + print '(A,I0,A,ES16.6,A,ES16.6)', ' moving_ibm = ', patch_ib(near1_id)%moving_ibm, ' mass=', & + & patch_ib(near1_id)%mass, ' moment=', patch_ib(near1_id)%moment + end if + if (near2_id > 0) then + print '(A,I0,A,ES16.6,A,ES16.6,A,3(ES16.6,1X))', ' 2nd nearest particle id=', near2_id, ' dist=', near2_dist, & + & ' gap=', near2_dist - patch_ib(near2_id)%radius, ' vel=', patch_ib(near2_id)%vel + print '(A,3(ES16.6,1X))', ' centroid = ', patch_ib(near2_id)%x_centroid, patch_ib(near2_id)%y_centroid, & + & patch_ib(near2_id)%z_centroid + print '(A,3(ES16.6,1X))', ' angular_vel = ', patch_ib(near2_id)%angular_vel + print '(A,3(ES16.6,1X))', ' force = ', patch_ib(near2_id)%force + end if + end if + + ! Dump a small x-neighborhood around the violating cell (reaching into the ghost/halo region on either side) to + ! distinguish a sharp discontinuity at a processor boundary - the signature of stale or corrupted halo/IB state - + ! from a smoothly diverging field, which indicates a genuine physical/numerical instability. + print '(A)', ' x-neighborhood (dj, rho, pres, vel) around violating cell:' + do j = max(-buff_size, j_hit - 3), min(m + buff_size, j_hit + 3) + call s_compute_cell_state(q_prim_vf, pres, rho, gamma, pi_inf, Re, alpha, alpha_rho, vel, vel_sum, qv, j, k_hit, l_hit) + print '(A,I0,A,ES16.6,A,ES16.6,A,3(ES16.6,1X))', ' dj=', j - j_hit, ' rho=', rho, ' pres=', pres, ' vel=', vel + end do + + call flush (6) + + end subroutine s_report_icfl_violation + !> Write grid and conservative variable data files in serial format impure subroutine s_write_serial_data_files(q_cons_vf, q_T_sf, q_prim_vf, t_step, bc_type, beta) diff --git a/src/simulation/m_derived_variables.fpp b/src/simulation/m_derived_variables.fpp index ec0def1d36..bb0b25738a 100644 --- a/src/simulation/m_derived_variables.fpp +++ b/src/simulation/m_derived_variables.fpp @@ -37,14 +37,17 @@ contains ! higher than fourth-order accuracy coefficients are wanted, the formulae required to compute these coefficients will have ! to be implemented in the subroutine s_compute_finite_difference_coefficients. - ! Allocating centered finite-difference coefficients + ! Allocating centered finite-difference coefficients. The coefficient (second) index is extended by fd_number beyond + ! the interior on each side: s_compute_ib_forces evaluates the viscous-stress stencil centered on ghost-adjacent + ! cells (i+l for l in -fd_number:fd_number) when an IB sits near a domain boundary, so the coefficient array must + ! cover those centers too, not just the interior 0:m. if (probe_wrt .or. ib) then - @:ALLOCATE(fd_coeff_x(-fd_number:fd_number, 0:m)) + @:ALLOCATE(fd_coeff_x(-fd_number:fd_number,-fd_number:m + fd_number)) if (n > 0) then - @:ALLOCATE(fd_coeff_y(-fd_number:fd_number, 0:n)) + @:ALLOCATE(fd_coeff_y(-fd_number:fd_number,-fd_number:n + fd_number)) end if if (p > 0) then - @:ALLOCATE(fd_coeff_z(-fd_number:fd_number, 0:p)) + @:ALLOCATE(fd_coeff_z(-fd_number:fd_number,-fd_number:p + fd_number)) end if @:ALLOCATE(accel_mag(0:m, 0:n, 0:p)) @@ -67,7 +70,8 @@ contains if (proc_rank == 0 .and. probe_wrt) then call s_open_probe_files() end if - ! Computing centered finite difference coefficients + ! Computing centered finite difference coefficients (s_compute_finite_difference_coefficients always extends + ! fd_number beyond the interior on each side; the allocation above matches) call s_compute_finite_difference_coefficients(m, x_cc, fd_coeff_x, buff_size, fd_number, fd_order) $:GPU_UPDATE(device='[fd_coeff_x]') diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index 9f315c5cec..e72afb668c 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -675,6 +675,7 @@ contains particle_cloud(i)%moving_ibm = 0 particle_cloud(i)%seed = 0 particle_cloud(i)%cloud_geometry = 1 + particle_cloud(i)%shell_axis = 3 particle_cloud(i)%packing_method = dflt_int particle_cloud(i)%periodic = 0 end do diff --git a/src/simulation/m_hypoelastic.fpp b/src/simulation/m_hypoelastic.fpp index 187fa70ce4..9da6817e8c 100644 --- a/src/simulation/m_hypoelastic.fpp +++ b/src/simulation/m_hypoelastic.fpp @@ -59,12 +59,13 @@ contains end do $:GPU_UPDATE(device='[Gs_hypo]') - @:ALLOCATE(fd_coeff_x_hypo(-fd_number:fd_number, 0:m)) + ! s_compute_finite_difference_coefficients always extends fd_number beyond the interior on each side + @:ALLOCATE(fd_coeff_x_hypo(-fd_number:fd_number,-fd_number:m + fd_number)) if (n > 0) then - @:ALLOCATE(fd_coeff_y_hypo(-fd_number:fd_number, 0:n)) + @:ALLOCATE(fd_coeff_y_hypo(-fd_number:fd_number,-fd_number:n + fd_number)) end if if (p > 0) then - @:ALLOCATE(fd_coeff_z_hypo(-fd_number:fd_number, 0:p)) + @:ALLOCATE(fd_coeff_z_hypo(-fd_number:fd_number,-fd_number:p + fd_number)) end if ! Computing centered finite difference coefficients diff --git a/src/simulation/m_ib_patches.fpp b/src/simulation/m_ib_patches.fpp index 0a845c8fe1..84e3df5782 100644 --- a/src/simulation/m_ib_patches.fpp +++ b/src/simulation/m_ib_patches.fpp @@ -228,38 +228,53 @@ contains ! rotate the frame into the IB's coordinates xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local) - ! perform the interior check for the patch geometry of this IB + ! perform the interior check for the patch geometry of this IB. Writes to ib_markers use + ! an atomic max (not a plain assignment) because this loop is parallel over patch_id: the + ! soft-sphere collision model allows particles to physically interpenetrate by design, so + ! two different patches can both claim the same overlapping cell here. A plain write would + ! be a data race with a nondeterministic winner; the atomic max makes the higher + ! encoded_patch_id win consistently every time, regardless of thread scheduling. if (patch_ib(patch_id)%geometry == 8) then ! sphere geometry radius = patch_ib(patch_id)%radius - if (f_is_inside_sphere(xyz_local(1), xyz_local(2), xyz_local(3), & - & radius)) ib_markers%sf(i, j, k) = encoded_patch_id + if (f_is_inside_sphere(xyz_local(1), xyz_local(2), xyz_local(3), radius)) then + $:GPU_ATOMIC(atomic='update') + ib_markers%sf(i, j, k) = max(ib_markers%sf(i, j, k), encoded_patch_id) + end if else if (patch_ib(patch_id)%geometry == 9) then ! cuboid geometry length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, & & patch_ib(patch_id)%length_z] - if (f_is_inside_cuboid(xyz_local(1), xyz_local(2), xyz_local(3), & - & length)) ib_markers%sf(i, j, k) = encoded_patch_id + if (f_is_inside_cuboid(xyz_local(1), xyz_local(2), xyz_local(3), length)) then + $:GPU_ATOMIC(atomic='update') + ib_markers%sf(i, j, k) = max(ib_markers%sf(i, j, k), encoded_patch_id) + end if else if (patch_ib(patch_id)%geometry == 10) then ! cylinder geometry radius = patch_ib(patch_id)%radius if (f_is_inside_cylinder(xyz_local(2), xyz_local(3), xyz_local(1), radius, & - & patch_ib(patch_id)%length_x)) ib_markers%sf(i, j, k) = encoded_patch_id + & patch_ib(patch_id)%length_x)) then + $:GPU_ATOMIC(atomic='update') + ib_markers%sf(i, j, k) = max(ib_markers%sf(i, j, k), encoded_patch_id) + end if else if (patch_ib(patch_id)%geometry == 11) then ! 3D airfoil geometry airfoil_id = patch_ib(patch_id)%airfoil_id xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset if (f_is_inside_airfoil(xyz_local(1), xyz_local(2), xyz_local(3), & - & patch_ib(patch_id)%length_z, airfoil_id)) ib_markers%sf(i, j, & - & k) = encoded_patch_id + & patch_ib(patch_id)%length_z, airfoil_id)) then + $:GPU_ATOMIC(atomic='update') + ib_markers%sf(i, j, k) = max(ib_markers%sf(i, j, k), encoded_patch_id) + end if else if (patch_ib(patch_id)%geometry == 12) then ! STL model geometry xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset model_id = patch_ib(patch_id)%model_id eta = f_model_is_inside(gpu_ntrs(model_id), model_id, xyz_local) if (eta > stl_models(model_id)%model_threshold) then - ib_markers%sf(i, j, k) = encoded_patch_id + $:GPU_ATOMIC(atomic='update') + ib_markers%sf(i, j, k) = max(ib_markers%sf(i, j, k), encoded_patch_id) end if end if end do @@ -296,36 +311,50 @@ contains ! rotate the frame into the IB's coordinates xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local) - ! perform the interior check for the patch geometry of this IB + ! perform the interior check for the patch geometry of this IB. Writes to ib_markers use an + ! atomic max (not a plain assignment) because this loop is parallel over patch_id: the + ! soft-sphere collision model allows particles to physically interpenetrate by design, so two + ! different patches can both claim the same overlapping cell here. A plain write would be a data + ! race with a nondeterministic winner; the atomic max makes the higher encoded_patch_id win + ! consistently every time, regardless of thread scheduling. if (patch_ib(patch_id)%geometry == 2) then ! circular geometries radius = patch_ib(patch_id)%radius - if (f_is_inside_cylinder(xyz_local(1), xyz_local(2), 0._wp, radius, 0._wp)) ib_markers%sf(i, & - & j, 0) = encoded_patch_id + if (f_is_inside_cylinder(xyz_local(1), xyz_local(2), 0._wp, radius, 0._wp)) then + $:GPU_ATOMIC(atomic='update') + ib_markers%sf(i, j, 0) = max(ib_markers%sf(i, j, 0), encoded_patch_id) + end if else if (patch_ib(patch_id)%geometry == 3) then ! rectangular geometries length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, 0._wp] - if (f_is_inside_cuboid(xyz_local(1), xyz_local(2), xyz_local(3), length)) ib_markers%sf(i, j, & - & 0) = encoded_patch_id + if (f_is_inside_cuboid(xyz_local(1), xyz_local(2), xyz_local(3), length)) then + $:GPU_ATOMIC(atomic='update') + ib_markers%sf(i, j, 0) = max(ib_markers%sf(i, j, 0), encoded_patch_id) + end if else if (patch_ib(patch_id)%geometry == 4) then ! 2D airfoil geometry airfoil_id = patch_ib(patch_id)%airfoil_id xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset - if (f_is_inside_airfoil(xyz_local(1), xyz_local(2), 0._wp, 0._wp, & - & airfoil_id)) ib_markers%sf(i, j, 0) = encoded_patch_id + if (f_is_inside_airfoil(xyz_local(1), xyz_local(2), 0._wp, 0._wp, airfoil_id)) then + $:GPU_ATOMIC(atomic='update') + ib_markers%sf(i, j, 0) = max(ib_markers%sf(i, j, 0), encoded_patch_id) + end if else if (patch_ib(patch_id)%geometry == 5) then ! STL model geometry xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset model_id = patch_ib(patch_id)%model_id eta = f_model_is_inside(gpu_ntrs(model_id), model_id, xyz_local) if (eta > stl_models(model_id)%model_threshold) then - ib_markers%sf(i, j, 0) = encoded_patch_id + $:GPU_ATOMIC(atomic='update') + ib_markers%sf(i, j, 0) = max(ib_markers%sf(i, j, 0), encoded_patch_id) end if else if (patch_ib(patch_id)%geometry == 6) then ! ellipse geometry length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, 0._wp] - if (f_is_inside_ellipse(xyz_local(1), xyz_local(2), length)) ib_markers%sf(i, j, & - & 0) = encoded_patch_id + if (f_is_inside_ellipse(xyz_local(1), xyz_local(2), length)) then + $:GPU_ATOMIC(atomic='update') + ib_markers%sf(i, j, 0) = max(ib_markers%sf(i, j, 0), encoded_patch_id) + end if end if end do end do @@ -572,10 +601,12 @@ contains end if ! completely skip patches whose bounding box does not overlap this rank's domain - outside_domain = bbox_min(1) > x_cc(m + gp_layers + 1) .or. bbox_max(1) < x_cc(-gp_layers - 1) .or. bbox_min(2) > y_cc(n & - & + gp_layers + 1) .or. bbox_max(2) < y_cc(-gp_layers - 1) + ! Markers cover the full halo: image-point stencils of ghost points near a rank boundary read markers up to + ! 2*gp_layers+1 cells into it, and an undrawn cell there reads as fluid. + outside_domain = bbox_min(1) > x_cc(m + buff_size) .or. bbox_max(1) < x_cc(-buff_size) .or. bbox_min(2) > y_cc(n & + & + buff_size) .or. bbox_max(2) < y_cc(-buff_size) if (num_dims == 3) then - outside_domain = outside_domain .or. bbox_min(3) > z_cc(p + gp_layers + 1) .or. bbox_max(3) < z_cc(-gp_layers - 1) + outside_domain = outside_domain .or. bbox_min(3) > z_cc(p + buff_size) .or. bbox_max(3) < z_cc(-buff_size) end if if (outside_domain) then @@ -585,12 +616,12 @@ contains return end if - il = -gp_layers - 1 - jl = -gp_layers - 1 - kl = -gp_layers - 1 - ir = m + gp_layers + 1 - jr = n + gp_layers + 1 - kr = p + gp_layers + 1 + il = -buff_size + jl = -buff_size + kl = -buff_size + ir = m + buff_size + jr = n + buff_size + kr = p + buff_size call get_indices_from_bounds(bbox_min(1), bbox_max(1), x_cc, il, ir) call get_indices_from_bounds(bbox_min(2), bbox_max(2), y_cc, jl, jr) if (num_dims == 3) call get_indices_from_bounds(bbox_min(3), bbox_max(3), z_cc, kl, kr) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 543a0fa6a6..665dbb5de0 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -33,6 +33,11 @@ module m_ibm type(integer_field), public :: ib_markers $:GPU_DECLARE(create='[ib_markers]') + !> 1 at ghost points whose image point reaches the fluid, 0 elsewhere. Ghost points with a buried image point average over these + !! neighbors, and never over each other. + type(integer_field) :: corrected_gps + $:GPU_DECLARE(create='[corrected_gps]') + type(ghost_point), dimension(:), allocatable :: ghost_points $:GPU_DECLARE(create='[ghost_points]') @@ -56,11 +61,14 @@ contains if (p > 0) then @:ALLOCATE(ib_markers%sf(-buff_size:m+buff_size, -buff_size:n+buff_size, -buff_size:p+buff_size)) + @:ALLOCATE(corrected_gps%sf(-buff_size:m+buff_size, -buff_size:n+buff_size, -buff_size:p+buff_size)) else @:ALLOCATE(ib_markers%sf(-buff_size:m+buff_size, -buff_size:n+buff_size, 0:0)) + @:ALLOCATE(corrected_gps%sf(-buff_size:m+buff_size, -buff_size:n+buff_size, 0:0)) end if @:ACC_SETUP_SFs(ib_markers) + @:ACC_SETUP_SFs(corrected_gps) $:GPU_ENTER_DATA(copyin='[num_gps]') @@ -112,7 +120,8 @@ contains ! recompute the new ib_patch locations ib_markers%sf = 0._wp - $:GPU_UPDATE(device='[ib_markers%sf]') + corrected_gps%sf = 0 + $:GPU_UPDATE(device='[ib_markers%sf, corrected_gps%sf]') call s_apply_ib_patches(ib_markers) $:GPU_UPDATE(host='[ib_markers%sf]') do i = 1, num_ibs @@ -145,28 +154,21 @@ contains end subroutine s_ibm_setup - subroutine s_compute_ghost_point_pressure(gp, gp_patch_id, alpha_rho_IP, pres_IP, pres_GP) + !> Pressure correction for a moving IB, accounting for the acceleration of the boundary surface. Clamped both ways: the + !! linearization it comes from holds only while the correction is order one, and an unbounded one drives the ghost state to + !! vacuum. + subroutine s_compute_ghost_point_pressure(gp, gp_patch_id, rho, pres_IP, pres_GP) $:GPU_ROUTINE(parallelism='[seq]') type(ghost_point), intent(in) :: gp integer, intent(in) :: gp_patch_id - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in) :: alpha_rho_IP - #:else - real(wp), dimension(num_fluids), intent(in) :: alpha_rho_IP - #:endif - real(wp), intent(in) :: pres_IP - real(wp), intent(out) :: pres_GP - integer :: q !< Iterator variable + real(wp), intent(in) :: rho, pres_IP + real(wp), intent(out) :: pres_GP - pres_GP = 0._wp - $:GPU_LOOP(parallelism='[seq]') - do q = 1, num_fluids - ! Pressure correction for moving IB: accounts for acceleration of IB surface - pres_GP = pres_GP + pres_IP/(1._wp - 2._wp*abs(gp%levelset*alpha_rho_IP(q)/pres_IP) & - & *dot_product(patch_ib(gp_patch_id)%force/patch_ib(gp_patch_id)%mass, gp%levelset_norm)) - end do + pres_GP = pres_IP/min(max(1._wp - 2._wp*abs(gp%levelset) & + & *rho/pres_IP*dot_product(patch_ib(gp_patch_id)%force/patch_ib(gp_patch_id)%mass, & + & gp%levelset_norm), 5.e-1_wp), 2._wp) end subroutine s_compute_ghost_point_pressure @@ -237,6 +239,8 @@ contains type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf !< Primitive Variables real(stp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), optional, intent(inout) :: pb_in, mv_in integer :: i, j, k, l, q, r !< Iterator variables + integer :: jj, kk, ll !< Neighbor iterators + integer :: rad, rad_z, num_nbrs !< Neighbor stencil radius and population integer :: patch_id, patch_id_temp !< Patch ID of ghost point real(wp) :: rho, gamma, pi_inf, dyn_pres !< Mixture variables real(wp) :: vel_sum_g, E_ghost !< Ghost-point velocity magnitude and energy @@ -249,14 +253,14 @@ contains #:if not MFC_CASE_OPTIMIZATION and USING_AMD real(wp), dimension(3) :: Gs - real(wp), dimension(3) :: alpha_rho_IP, alpha_IP + real(wp), dimension(3) :: alpha_rho_IP, alpha_IP, alpha_rho_GP real(wp), dimension(3) :: r_IP, v_IP, pb_IP, mv_IP real(wp), dimension(18) :: nmom_IP real(wp), dimension(12) :: presb_IP, massv_IP real(wp), dimension(${AMD_NUM_SPECIES_MAX}$) :: Ys_IP #:else real(wp), dimension(num_fluids) :: Gs - real(wp), dimension(num_fluids) :: alpha_rho_IP, alpha_IP + real(wp), dimension(num_fluids) :: alpha_rho_IP, alpha_IP, alpha_rho_GP real(wp), dimension(nb) :: r_IP, v_IP, pb_IP, mv_IP real(wp), dimension(nb*nmom) :: nmom_IP real(wp), dimension(nb*nnode) :: presb_IP, massv_IP @@ -270,10 +274,11 @@ contains real(wp), dimension(3) :: vel_g !< Velocity of GP real(wp), dimension(3) :: radial_vector !< vector from centroid to ghost point real(wp) :: nbub + real(wp) :: buf, buf_prim type(ghost_point) :: gp ! set the Moving IBM interior conservative variables - $:GPU_PARALLEL_LOOP(private='[i, j, k, patch_id, rho]', collapse=3) + $:GPU_PARALLEL_LOOP(private='[i, j, k, patch_id, rho, patch_id_temp]', collapse=3) do l = 0, p do k = 0, n do j = 0, m @@ -282,12 +287,7 @@ contains call s_decode_patch_periodicity(patch_id, patch_id_temp) call s_get_neighborhood_idx(patch_id_temp, patch_id) if (patch_id > 0) then - ! Placeholder low pressure inside the IB solid. Skip it with - ! chemistry on: it would force an unphysical temperature - ! (P=1 Pa at the ambient density -> T~0.01 K), which the - ! Cantera temperature/transport evaluation (run grid-wide - ! before the IB mask is applied) cannot handle -> NaN/hang. - ! The interior is masked from the RHS regardless. + ! skip pressure correction with chemistry to prevent unphysical pressure if (.not. chemistry) q_prim_vf(eqn_idx%E)%sf(j, k, l) = 1._wp rho = 0._wp do i = 1, num_fluids @@ -307,12 +307,13 @@ contains $:END_GPU_PARALLEL_LOOP() if (num_gps > 0) then - $:GPU_PARALLEL_LOOP(private='[i, physical_loc, dyn_pres, alpha_rho_IP, alpha_IP, pres_IP, pres_GP, vel_IP, vel_g, & - & r_IP, v_IP, pb_IP, mv_IP, nmom_IP, presb_IP, massv_IP, rho, gamma, pi_inf, Re_K, G_K, Gs, gp, & - & radial_vector, j, k, l, q, qv_K, c_IP, nbub, patch_id, Ys_IP, T_IP, mw_IP, e_IP, vel_sum_g, & - & E_ghost, alpha_q, alpha_rho_q, e_q]') + $:GPU_PARALLEL_LOOP(private='[i, physical_loc, dyn_pres, alpha_rho_IP, alpha_IP, alpha_rho_GP, pres_IP, pres_GP, & + & vel_IP, vel_g, r_IP, v_IP, pb_IP, mv_IP, nmom_IP, presb_IP, massv_IP, rho, gamma, pi_inf, Re_K, & + & G_K, Gs, gp, radial_vector, j, k, l, q, qv_K, c_IP, nbub, patch_id, Ys_IP, T_IP, mw_IP, e_IP, & + & vel_sum_g, E_ghost, alpha_q, alpha_rho_q, e_q]', present='[ghost_points]') do i = 1, num_gps gp = ghost_points(i) + if (.not. gp%interp_valid) cycle j = gp%loc(1) k = gp%loc(2) l = gp%loc(3) @@ -351,27 +352,6 @@ contains alpha_rho_IP(1) = pres_IP*mw_IP/(T_IP*gas_constant) end if - dyn_pres = 0._wp - - ! Set q_prim_vf params at GP so that mixture vars calculated properly - $:GPU_LOOP(parallelism='[seq]') - do q = 1, num_fluids - q_prim_vf(q)%sf(j, k, l) = alpha_rho_IP(q) - q_prim_vf(eqn_idx%adv%beg + q - 1)%sf(j, k, l) = alpha_IP(q) - end do - - if (surface_tension) then - q_prim_vf(eqn_idx%c)%sf(j, k, l) = c_IP - end if - - ! set the pressure - if (patch_ib(patch_id)%moving_ibm <= 1) then - q_prim_vf(eqn_idx%E)%sf(j, k, l) = pres_IP - else - call s_compute_ghost_point_pressure(gp, patch_id, alpha_rho_IP, pres_IP, pres_GP) - q_prim_vf(eqn_idx%E)%sf(j, k, l) = pres_GP - end if - ! If in simulation, use acc mixture subroutines if (hypoelasticity) then call s_convert_species_to_mixture_variables_kernel(rho, gamma, pi_inf, qv_K, alpha_IP, alpha_rho_IP, Re_K, & @@ -380,6 +360,34 @@ contains call s_convert_species_to_mixture_variables_kernel(rho, gamma, pi_inf, qv_K, alpha_IP, alpha_rho_IP, Re_K) end if + if (surface_tension) q_prim_vf(eqn_idx%c)%sf(j, k, l) = c_IP + + ! set the pressure and density + if (patch_ib(patch_id)%moving_ibm <= 1) then + pres_GP = pres_IP + alpha_rho_GP = alpha_rho_IP + else + call s_compute_ghost_point_pressure(gp, patch_id, rho, pres_IP, pres_GP) + + ! The adiabatic wall condition T_GP = T_IP the correction is derived from also + ! fixes the ghost density: p + B = (n - 1)*cv*rho*T at both points under the one + ! temperature leaves each partial density carrying the pressure ratio. The volume + ! fractions are untouched, so only rho and qv move with it. + $:GPU_LOOP(parallelism='[seq]') + do q = 1, num_fluids + alpha_rho_GP(q) = alpha_rho_IP(q)*(pres_GP + isentrope_B(q))/(pres_IP + isentrope_B(q)) + end do + call s_compute_mixture_coefficients(alpha_rho_GP, alpha_IP, rho, gamma, pi_inf, qv_K) + end if + q_prim_vf(eqn_idx%E)%sf(j, k, l) = pres_GP + + ! Set q_prim_vf params at GP + $:GPU_LOOP(parallelism='[seq]') + do q = 1, num_fluids + q_prim_vf(q)%sf(j, k, l) = alpha_rho_GP(q) + q_prim_vf(eqn_idx%adv%beg + q - 1)%sf(j, k, l) = alpha_IP(q) + end do + ! get the vector that points from the centroid to the ghost radial_vector(1) = physical_loc(1) - (patch_ib(patch_id)%x_centroid + real(ghost_points(i)%x_periodicity, & & wp)*(glb_bounds(1)%end - glb_bounds(1)%beg)) @@ -397,6 +405,7 @@ contains $:GPU_LOOP(parallelism='[seq]') do q = eqn_idx%mom%beg, eqn_idx%mom%end q_cons_vf(q)%sf(j, k, l) = rho*vel_g(q - eqn_idx%mom%beg + 1) + q_prim_vf(q)%sf(j, k, l) = vel_g(q - eqn_idx%mom%beg + 1) vel_sum_g = vel_sum_g + vel_g(q - eqn_idx%mom%beg + 1)**2._wp end do dyn_pres = 5.e-1_wp*rho*vel_sum_g @@ -404,7 +413,7 @@ contains ! Set continuity and adv vars $:GPU_LOOP(parallelism='[seq]') do q = 1, num_fluids - q_cons_vf(q)%sf(j, k, l) = alpha_rho_IP(q) + q_cons_vf(q)%sf(j, k, l) = alpha_rho_GP(q) q_cons_vf(eqn_idx%adv%beg + q - 1)%sf(j, k, l) = alpha_IP(q) end do @@ -427,7 +436,7 @@ contains end do q_cons_vf(eqn_idx%E)%sf(j, k, l) = rho*e_IP + dyn_pres else - call s_compute_energy(pres_IP, alpha_rho_IP, alpha_IP, vel_sum_g, E_ghost) + call s_compute_energy(pres_GP, alpha_rho_GP, alpha_IP, vel_sum_g, E_ghost) q_cons_vf(eqn_idx%E)%sf(j, k, l) = E_ghost end if ! Set bubble vars @@ -474,13 +483,86 @@ contains $:GPU_LOOP(parallelism='[seq]') do q = eqn_idx%int_en%beg, eqn_idx%int_en%end alpha_q = alpha_IP(q - eqn_idx%int_en%beg + 1) - alpha_rho_q = alpha_rho_IP(q - eqn_idx%int_en%beg + 1) - call s_phase_internal_energy(pres_IP, alpha_q, alpha_rho_q, q - eqn_idx%int_en%beg + 1, e_q) + alpha_rho_q = alpha_rho_GP(q - eqn_idx%int_en%beg + 1) + call s_phase_internal_energy(pres_GP, alpha_q, alpha_rho_q, q - eqn_idx%int_en%beg + 1, e_q) q_cons_vf(q)%sf(j, k, l) = e_q end do end if end do $:END_GPU_PARALLEL_LOOP() + + ! A ghost point whose image point is buried in a neighboring IB has no fluid to mirror, so + ! it takes the average of the ghost points corrected above, growing the stencil until it + ! reaches one. Those neighbors already carry the wall condition, so the average does too. + $:GPU_PARALLEL_LOOP(private='[i, j, k, l, q, r, jj, kk, ll, gp, rad, rad_z, num_nbrs, buf, buf_prim]', & + & present='[ghost_points]') + do i = 1, num_gps + gp = ghost_points(i) + if (gp%interp_valid) cycle + j = gp%loc(1) + k = gp%loc(2) + l = gp%loc(3) + + num_nbrs = 0 + rad = 0 + rad_z = 0 + do while (num_nbrs == 0 .and. rad < gp_layers) + rad = rad + 1 + if (p /= 0) rad_z = rad + do jj = j - rad, j + rad + do kk = k - rad, k + rad + do ll = l - rad_z, l + rad_z + num_nbrs = num_nbrs + corrected_gps%sf(jj, kk, ll) + end do + end do + end do + end do + + ! Fully enclosed: nothing to average from, so leave the interior placeholder in place + if (num_nbrs == 0) cycle + + $:GPU_LOOP(parallelism='[seq]') + do q = 1, sys_size + buf = 0._wp + buf_prim = 0._wp + do jj = j - rad, j + rad + do kk = k - rad, k + rad + do ll = l - rad_z, l + rad_z + if (corrected_gps%sf(jj, kk, ll) == 1) then + buf = buf + q_cons_vf(q)%sf(jj, kk, ll) + buf_prim = buf_prim + q_prim_vf(q)%sf(jj, kk, ll) + end if + end do + end do + end do + q_cons_vf(q)%sf(j, k, l) = buf/real(num_nbrs, wp) + q_prim_vf(q)%sf(j, k, l) = buf_prim/real(num_nbrs, wp) + end do + + if (qbmm .and. .not. polytropic) then + $:GPU_LOOP(parallelism='[seq]') + do q = 1, nb + $:GPU_LOOP(parallelism='[seq]') + do r = 1, nnode + buf = 0._wp + buf_prim = 0._wp + do jj = j - rad, j + rad + do kk = k - rad, k + rad + do ll = l - rad_z, l + rad_z + if (corrected_gps%sf(jj, kk, ll) == 1) then + buf = buf + pb_in(jj, kk, ll, r, q) + buf_prim = buf_prim + mv_in(jj, kk, ll, r, q) + end if + end do + end do + end do + pb_in(j, k, l, r, q) = buf/real(num_nbrs, wp) + mv_in(j, k, l, r, q) = buf_prim/real(num_nbrs, wp) + end do + end do + end if + end do + $:END_GPU_PARALLEL_LOOP() end if end subroutine s_ibm_correct_state @@ -734,10 +816,9 @@ contains real(wp), dimension(2, 2, 2) :: eta type(ghost_point) :: gp integer :: q, i, j, k, ii, jj, kk !< Grid indexes and iterators - integer :: patch_id logical :: is_cell_center - $:GPU_PARALLEL_LOOP(private='[q, i, j, k, ii, jj, kk, dist, buf, gp, interp_coeffs, eta, alpha, patch_id, is_cell_center]') + $:GPU_PARALLEL_LOOP(private='[q, i, j, k, ii, jj, kk, dist, buf, gp, interp_coeffs, eta, alpha, is_cell_center]') do q = 1, num_gps gp = ghost_points_in(q) ! Get the interpolation points @@ -771,14 +852,16 @@ contains is_cell_center = .false. check_is_cell_center: do ii = 0, 1 do jj = 0, 1 + ! A coincident solid cell leaves the point invalid: falling through to the eta + ! branch with dist = 0 would give alpha*eta = 0*Inf = NaN. if (dist(ii + 1, jj + 1, 1) <= 1.e-16_wp) then - interp_coeffs(ii + 1, jj + 1, 1) = 1._wp + if (ib_markers%sf(i + ii, j + jj, k) == 0) interp_coeffs(ii + 1, jj + 1, 1) = 1._wp is_cell_center = .true. exit check_is_cell_center else if (p /= 0) then if (dist(ii + 1, jj + 1, 2) <= 1.e-16_wp) then - interp_coeffs(ii + 1, jj + 1, 2) = 1._wp + if (ib_markers%sf(i + ii, j + jj, k + 1) == 0) interp_coeffs(ii + 1, jj + 1, 2) = 1._wp is_cell_center = .true. exit check_is_cell_center end if @@ -790,7 +873,6 @@ contains if (.not. is_cell_center) then ! if we are not arbitrarily close, interpolate alpha = 1._wp - patch_id = gp%ib_patch_id if (ib_markers%sf(i, j, k) /= 0) alpha(1, 1, 1) = 0._wp if (ib_markers%sf(i + 1, j, k) /= 0) alpha(2, 1, 1) = 0._wp if (ib_markers%sf(i, j + 1, k) /= 0) alpha(1, 2, 1) = 0._wp @@ -799,12 +881,7 @@ contains if (p == 0) then eta(:,:,1) = 1._wp/dist(:,:,1)**2 buf = sum(alpha(:,:,1)*eta(:,:,1)) - if (buf > 0._wp) then - interp_coeffs(:,:,1) = alpha(:,:,1)*eta(:,:,1)/buf - else - buf = sum(eta(:,:,1)) - interp_coeffs(:,:,1) = eta(:,:,1)/buf - end if + if (buf > 0._wp) interp_coeffs(:,:,1) = alpha(:,:,1)*eta(:,:,1)/buf else if (ib_markers%sf(i, j, k + 1) /= 0) alpha(1, 1, 2) = 0._wp if (ib_markers%sf(i + 1, j, k + 1) /= 0) alpha(2, 1, 2) = 0._wp @@ -812,17 +889,15 @@ contains if (ib_markers%sf(i + 1, j + 1, k + 1) /= 0) alpha(2, 2, 2) = 0._wp eta = 1._wp/dist**2 buf = sum(alpha*eta) - - if (buf > 0._wp) then - interp_coeffs = alpha*eta/buf - else - buf = sum(eta) - interp_coeffs = eta/buf - end if + if (buf > 0._wp) interp_coeffs = alpha*eta/buf end if end if + ! An image point buried in a neighboring IB gets no weights at all: its ghost point is + ! averaged from its own neighbors instead, in s_ibm_correct_state. ghost_points_in(q)%interp_coeffs = interp_coeffs + ghost_points_in(q)%interp_valid = any(interp_coeffs > 0._wp) + if (ghost_points_in(q)%interp_valid) corrected_gps%sf(gp%loc(1), gp%loc(2), gp%loc(3)) = 1 end do $:END_GPU_PARALLEL_LOOP() @@ -960,15 +1035,16 @@ contains impure subroutine s_update_mib(num_ibs) integer, intent(in) :: num_ibs - integer :: i, j, k, z_gp_layers + integer :: i, j, k, z_buff_size call nvtxStartRange("UPDATE-MIBM") - ! Clears the existing immersed boundary indices - z_gp_layers = 0; if (p /= 0) z_gp_layers = gp_layers + 1 + ! Clears the existing immersed boundary indices over the same halo extent s_get_bounding_indices draws them + z_buff_size = 0; if (p /= 0) z_buff_size = buff_size $:GPU_PARALLEL_LOOP(private='[i, j, k]') - do i = -gp_layers - 1, m + gp_layers + 1; do j = -gp_layers - 1, n + gp_layers + 1; do k = -z_gp_layers, p + z_gp_layers + do i = -buff_size, m + buff_size; do j = -buff_size, n + buff_size; do k = -z_buff_size, p + z_buff_size ib_markers%sf(i, j, k) = 0._wp + corrected_gps%sf(i, j, k) = 0 end do; end do; end do $:END_GPU_PARALLEL_LOOP() @@ -1222,10 +1298,12 @@ contains end do ! apply the summed forces - $:GPU_PARALLEL_LOOP(private='[i]', copyin='[forces, torques]') + $:GPU_PARALLEL_LOOP(private='[i, l]', copyin='[forces, torques]') do i = 1, num_ibs - patch_ib(i)%force(:) = forces(i,:) - patch_ib(i)%torque(:) = torques(i,:) + do l = 1, 3 + patch_ib(i)%force(l) = forces(i, l) + patch_ib(i)%torque(l) = torques(i, l) + end do end do $:END_GPU_PARALLEL_LOOP() @@ -1409,7 +1487,7 @@ contains 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 :: i, j, k, l, 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(:) @@ -1432,11 +1510,13 @@ contains do k = 1, min(2*ib_neighborhood_radius, num_procs_${X}$ - 1) ! send forces to +${X}$ neighbor; receive from -${X}$ neighbor. Add received values then pack_pos = 0 - $:GPU_PARALLEL_LOOP(private='[i]', copyin='[forces, torques]') + $:GPU_PARALLEL_LOOP(private='[i, l]', 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,:) + do l = 1, 3 + send_ft(l, i) = forces(i, l) + send_ft(l + 3, i) = torques(i, l) + end do end do $:END_GPU_PARALLEL_LOOP() $:GPU_UPDATE(host='[send_ids, send_ft]') @@ -1453,15 +1533,17 @@ contains & 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_UPDATE(device='[recv_ids(1:recv_count), recv_ft(:, 1:recv_count)]') - $:GPU_PARALLEL_LOOP(private='[i, j]', copy='[forces, torques]') + $:GPU_PARALLEL_LOOP(private='[i, j, l]', copy='[forces, torques]') 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) + do l = 1, 3 + forces(j, l) = forces(j, l) + recv_ft(l, i) - recv_forces_snap(j, l) + torques(j, l) = torques(j, l) + recv_ft(l + 3, i) - recv_torques_snap(j, l) + recv_forces_snap(j, l) = recv_ft(l, i) + recv_torques_snap(j, l) = recv_ft(l + 3, i) + end do end if end do $:END_GPU_PARALLEL_LOOP() @@ -1479,11 +1561,13 @@ contains do k = 1, min(2*ib_neighborhood_radius, num_procs_${X}$ - 1) pack_pos = 0 - $:GPU_PARALLEL_LOOP(private='[i]', copyin='[forces, torques]') + $:GPU_PARALLEL_LOOP(private='[i, l]', 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,:) + do l = 1, 3 + send_ft(l, i) = forces(i, l) + send_ft(l + 3, i) = torques(i, l) + end do end do $:END_GPU_PARALLEL_LOOP() $:GPU_UPDATE(host='[send_ids, send_ft]') @@ -1499,12 +1583,14 @@ contains & 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_UPDATE(device='[recv_ids(1:recv_count), recv_ft(:, 1:recv_count)]') - $:GPU_PARALLEL_LOOP(private='[i, j]', copy='[forces, torques]') + $:GPU_PARALLEL_LOOP(private='[i, j, l]', 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) + do l = 1, 3 + forces(j, l) = recv_ft(l, i) + torques(j, l) = recv_ft(l + 3, i) + end do end if end do $:END_GPU_PARALLEL_LOOP() @@ -1565,13 +1651,16 @@ contains ! check if in local domain if (f_local_rank_owns_location(centroid, glb_bounds)) then local_output_idx = local_output_idx + 1 + @:PROHIBIT(local_output_idx > num_local_ibs_max, & + & "Too many IBs on a single processor rank. Modify case file or increase limit of num_local_ibs_max to resolve.") 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]') + ! num_ibs shrinks here, so refresh it with patch_ib: s_update_ib_lookup scatters over it on the device + $:GPU_UPDATE(device='[patch_ib, num_ibs]') call s_update_ib_lookup() ! Broadcast newly-owned patches to all neighborhood neighbors @@ -1608,12 +1697,14 @@ contains ! 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 + do dz = merge(-ib_neighborhood_radius, 0, num_dims == 3), merge(ib_neighborhood_radius, 0, num_dims == 3) + do dy = -ib_neighborhood_radius, ib_neighborhood_radius + do dx = -ib_neighborhood_radius, ib_neighborhood_radius 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) + ! one tag per offset in the (2R+1)^3 neighbourhood: a radix-3 encoding collides once R > 1 + tag = 200 + ((dx + ib_neighborhood_radius)*(2*ib_neighborhood_radius + 1) + (dy + ib_neighborhood_radius)) & + & *(2*ib_neighborhood_radius + 1) + (dz + ib_neighborhood_radius) recv_neighbor = ib_neighbor_ranks(-dx, -dy, -dz) recv_neighbor_list(nbr_idx) = MPI_PROC_NULL if (recv_neighbor < 0) cycle @@ -1625,11 +1716,13 @@ contains 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 + do dz = merge(-ib_neighborhood_radius, 0, num_dims == 3), merge(ib_neighborhood_radius, 0, num_dims == 3) + do dy = -ib_neighborhood_radius, ib_neighborhood_radius + do dx = -ib_neighborhood_radius, ib_neighborhood_radius if (dx == 0 .and. dy == 0 .and. dz == 0) cycle - tag = 200 + (dx + 1)*9 + (dy + 1)*3 + (dz + 1) + ! one tag per offset in the (2R+1)^3 neighbourhood: a radix-3 encoding collides once R > 1 + tag = 200 + ((dx + ib_neighborhood_radius)*(2*ib_neighborhood_radius + 1) + (dy + ib_neighborhood_radius)) & + & *(2*ib_neighborhood_radius + 1) + (dz + ib_neighborhood_radius) send_neighbor = ib_neighbor_ranks(dx, dy, dz) if (send_neighbor < 0) cycle nreqs = nreqs + 1 @@ -1641,7 +1734,7 @@ contains call MPI_WAITALL(nreqs, requests, MPI_STATUSES_IGNORE, ierr) ! Unpack all received buffers - do nbr_idx = 1, merge(26, 8, num_dims == 3) + do nbr_idx = 1, ((2*ib_neighborhood_radius + 1)**num_dims) - 1 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) @@ -1658,7 +1751,7 @@ contains end do deallocate (send_buf, recv_bufs) - $:GPU_UPDATE(device='[patch_ib]') + $:GPU_UPDATE(device='[patch_ib, num_ibs]') call s_update_ib_lookup() end if #endif @@ -1700,6 +1793,7 @@ contains integer :: i @:DEALLOCATE(ib_markers%sf) + @:DEALLOCATE(corrected_gps%sf) @:DEALLOCATE(ib_gbl_idx_lookup) do i = 1, num_ib_airfoils_max if (allocated(ib_airfoil_grids(i)%upper)) then diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp index 132ba0d757..f5fdfd17e3 100644 --- a/src/simulation/m_mpi_proxy.fpp +++ b/src/simulation/m_mpi_proxy.fpp @@ -215,6 +215,7 @@ contains call MPI_BCAST(particle_cloud(i)%moving_ibm, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) call MPI_BCAST(particle_cloud(i)%seed, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) call MPI_BCAST(particle_cloud(i)%cloud_geometry, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(particle_cloud(i)%shell_axis, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) call MPI_BCAST(particle_cloud(i)%packing_method, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) call MPI_BCAST(particle_cloud(i)%periodic, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) end do diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 1c2d93afb1..6458a62201 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -591,7 +591,7 @@ contains end if if (dt < dt_floor .and. cfl_adap_dt .and. proc_rank == 0) then - print *, "Delta t = ", dt + print *, "Delta t = ", dt, " limited by ", dt_limiter call s_mpi_abort("Delta t has become too small") end if end if diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index c8a6953ce2..469b6577fa 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -679,6 +679,7 @@ contains real(wp), dimension(5) :: dt_candidates_loc !< Rank-local dt candidates (ICFL, VCFL, CCFL, TCFL, collision cap) real(wp), dimension(5) :: dt_candidates_glb !< Global dt candidates (ICFL, VCFL, CCFL, TCFL, collision cap) real(wp) :: dt_prev + logical :: is_fluid_cell !< Cell lies outside every immersed boundary integer :: j, k, l !< Generic loop iterators integer :: fl !< Fluid loop iterator @@ -693,39 +694,46 @@ contains tcfl_dt_local = huge(1.0_wp) coll_dt_local = huge(1.0_wp) $:GPU_PARALLEL_LOOP(collapse=3, private='[vel, alpha, alpha_rho, Re, rho, vel_sum, pres, gamma, pi_inf, c, qv, fl, & - & max_dt]', reduction='[[icfl_dt_local, vcfl_dt_local, ccfl_dt_local, tcfl_dt_local]]', reductionOp='[min]') + & max_dt, is_fluid_cell]', reduction='[[icfl_dt_local, vcfl_dt_local, ccfl_dt_local, & + & tcfl_dt_local]]', reductionOp='[min]') do l = 0, p do k = 0, n do j = 0, m - if (igr) then - call s_compute_cell_state(q_cons_ts(1)%vf, pres, rho, gamma, pi_inf, Re, alpha, alpha_rho, vel, vel_sum, & - & qv, j, k, l) - else - call s_compute_cell_state(q_prim_vf, pres, rho, gamma, pi_inf, Re, alpha, alpha_rho, vel, vel_sum, qv, j, & - & k, l) - end if + ! Cells inside an immersed boundary hold ghost-derived, non-physical state and must not set the global dt. + is_fluid_cell = .true. + if (ib) is_fluid_cell = (ib_markers%sf(j, k, l) == 0) + + if (is_fluid_cell) then + if (igr) then + call s_compute_cell_state(q_cons_ts(1)%vf, pres, rho, gamma, pi_inf, Re, alpha, alpha_rho, vel, & + & vel_sum, qv, j, k, l) + else + call s_compute_cell_state(q_prim_vf, pres, rho, gamma, pi_inf, Re, alpha, alpha_rho, vel, vel_sum, & + & qv, j, k, l) + end if + + ! Compute mixture sound speed + call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, alpha, c, alpha_rho) + + if (any_non_newtonian) then + Re(1) = 0._wp + do fl = 1, num_fluids + if (is_non_newtonian(fl)) then + Re(1) = Re(1) + alpha(fl)*hb_mu_max(fl) + else + Re(1) = Re(1) + alpha(fl)*fluid_inv_re(fl) + end if + end do + Re(1) = 1._wp/max(Re(1), sgm_eps) + end if - ! Compute mixture sound speed - call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, alpha, c, alpha_rho) + call s_compute_dt_from_cfl(vel, c, max_dt, rho, Re, alpha, alpha_rho, j, k, l) - if (any_non_newtonian) then - Re(1) = 0._wp - do fl = 1, num_fluids - if (is_non_newtonian(fl)) then - Re(1) = Re(1) + alpha(fl)*hb_mu_max(fl) - else - Re(1) = Re(1) + alpha(fl)*fluid_inv_re(fl) - end if - end do - Re(1) = 1._wp/max(Re(1), sgm_eps) + icfl_dt_local = min(icfl_dt_local, max_dt(1)) + vcfl_dt_local = min(vcfl_dt_local, max_dt(2)) + ccfl_dt_local = min(ccfl_dt_local, max_dt(3)) + tcfl_dt_local = min(tcfl_dt_local, max_dt(4)) end if - - call s_compute_dt_from_cfl(vel, c, max_dt, rho, Re, alpha, alpha_rho, j, k, l) - - icfl_dt_local = min(icfl_dt_local, max_dt(1)) - vcfl_dt_local = min(vcfl_dt_local, max_dt(2)) - ccfl_dt_local = min(ccfl_dt_local, max_dt(3)) - tcfl_dt_local = min(tcfl_dt_local, max_dt(4)) end do end do end do diff --git a/src/simulation/m_viscous.fpp b/src/simulation/m_viscous.fpp index 6084633f9f..37afd8988e 100644 --- a/src/simulation/m_viscous.fpp +++ b/src/simulation/m_viscous.fpp @@ -1127,33 +1127,28 @@ contains real(wp), dimension(1:3,1:3) :: velocity_gradient_tensor real(wp) :: divergence real(wp) :: mu_eff, gamma_dot_c - integer :: l, q !< iterators + integer :: l, q !< iterators integer :: fl integer :: r - integer :: i_fd, j_fd, k_fd !< sample clamped to the interior ! zero the viscous stress and collection of velocity derivatives viscous_stress_tensor = 0._wp velocity_gradient_tensor = 0._wp - ! fd_coeff_x/y/z are computed for interior cells only (0:m, 0:n, 0:p), but s_compute_ib_forces samples this routine - ! fd_number cells out from an interior cell, so a body near a domain boundary asks for a coefficient that was never - ! computed. Read the nearest interior cell's coefficients there: on a uniform grid they are the same, and on a - ! stretched one this is the stencil the boundary cell itself uses. - i_fd = min(max(i, 0), m) - j_fd = min(max(j, 0), n) - k_fd = min(max(k, 0), p) + ! s_compute_ib_forces centers this stencil up to fd_number cells outside the interior, so the coefficients are computed + ! that far beyond it too (s_compute_finite_difference_coefficients): every center read here has a real coefficient. + ! Clamping to the nearest interior cell instead would make a stretched-grid body's drag depend on the decomposition. ! compute the velocity gradient tensor with the same fd_order-respecting stencil as the stress-divergence outer derivative do l = 1, num_dims do r = -fd_number, fd_number velocity_gradient_tensor(l, 1) = velocity_gradient_tensor(l, 1) + fd_coeff_x(r, & - & i_fd)*q_prim_vf(eqn_idx%mom%beg + l - 1)%sf(i + r, j, k) + & i)*q_prim_vf(eqn_idx%mom%beg + l - 1)%sf(i + r, j, k) velocity_gradient_tensor(l, 2) = velocity_gradient_tensor(l, 2) + fd_coeff_y(r, & - & j_fd)*q_prim_vf(eqn_idx%mom%beg + l - 1)%sf(i, j + r, k) + & j)*q_prim_vf(eqn_idx%mom%beg + l - 1)%sf(i, j + r, k) if (num_dims == 3) then velocity_gradient_tensor(l, 3) = velocity_gradient_tensor(l, 3) + fd_coeff_z(r, & - & k_fd)*q_prim_vf(eqn_idx%mom%beg + l - 1)%sf(i, j, k + r) + & k)*q_prim_vf(eqn_idx%mom%beg + l - 1)%sf(i, j, k + r) end if end do end do diff --git a/tests/127A967A/golden-metadata.txt b/tests/127A967A/golden-metadata.txt index dcd99ff8f1..25bd38f250 100644 --- a/tests/127A967A/golden-metadata.txt +++ b/tests/127A967A/golden-metadata.txt @@ -1,19 +1,19 @@ -This file was created on 2025-12-04 18:03:33.497282. +This file was created on 2026-09-22 03:47:43.532521. mfc.sh: Invocation: test --only 127A967A --generate - Lock: mpi=Yes & gpu=No & debug=No & gcov=No & unified=No & single=No & mixed=No & fastmath=No - Git: cfe671402d160abd63c80750278f59293e426f85 on forces-via-volume-integrals (dirty) + Lock: mpi=Yes & gpu=No & debug=No & reldebug=No & gcov=No & unified=No & single=No & mixed=No & fastmath=No + Git: 414eba29941cedcf1c9bfd47781150ef7064ba4e on debug-ibm-stability (dirty) post_process: CMake Configuration: - CMake v4.0.3 on oppenheimer + CMake v3.28.3 on schwarzschild - C : AppleClang v17.0.0.17000013 (/usr/bin/cc) - Fortran : GNU v15.1.0 (/opt/homebrew/bin/gfortran) + C : NVHPC v23.11.0 (/opt/nvidia/hpc_sdk/Linux_x86_64/23.11/compilers/bin/nvc) + Fortran : NVHPC v23.11.0 (/opt/nvidia/hpc_sdk/Linux_x86_64/23.11/compilers/bin/nvfortran) PRE_PROCESS : OFF SIMULATION : OFF @@ -26,31 +26,31 @@ post_process: OpenACC : OFF OpenMP : OFF - Fypp : /Users/dan/Documents/repos/MFC/build/venv/bin/fypp + Fypp : /home/dan/Documents/repos/MFC/build/venv/bin/fypp Doxygen : Build Type : Release Configuration Environment: - CC : /usr/bin/cc - CXX : /usr/bin/c++ - FC : /opt/homebrew/bin/gfortran + CC : /opt/nvidia/hpc_sdk/Linux_x86_64/23.11/compilers/bin/nvc + CXX : /opt/nvidia/hpc_sdk/Linux_x86_64/23.11/compilers/bin/nvc++ + FC : /opt/nvidia/hpc_sdk/Linux_x86_64/23.11/compilers/bin/nvfortran OMPI_CC : OMPI_CXX : OMPI_FC : -pre_process: +simulation: CMake Configuration: - CMake v4.0.3 on oppenheimer + CMake v3.28.3 on schwarzschild - C : AppleClang v17.0.0.17000013 (/usr/bin/cc) - Fortran : GNU v15.1.0 (/opt/homebrew/bin/gfortran) + C : NVHPC v23.11.0 (/opt/nvidia/hpc_sdk/Linux_x86_64/23.11/compilers/bin/nvc) + Fortran : NVHPC v23.11.0 (/opt/nvidia/hpc_sdk/Linux_x86_64/23.11/compilers/bin/nvfortran) - PRE_PROCESS : ON - SIMULATION : OFF + PRE_PROCESS : OFF + SIMULATION : ON POST_PROCESS : OFF SYSCHECK : OFF DOCUMENTATION : OFF @@ -60,16 +60,16 @@ pre_process: OpenACC : OFF OpenMP : OFF - Fypp : /Users/dan/Documents/repos/MFC/build/venv/bin/fypp + Fypp : /home/dan/Documents/repos/MFC/build/venv/bin/fypp Doxygen : Build Type : Release Configuration Environment: - CC : /usr/bin/cc - CXX : /usr/bin/c++ - FC : /opt/homebrew/bin/gfortran + CC : nvc + CXX : nvc++ + FC : nvfortran OMPI_CC : OMPI_CXX : OMPI_FC : @@ -78,10 +78,10 @@ syscheck: CMake Configuration: - CMake v4.0.3 on oppenheimer + CMake v3.28.3 on schwarzschild - C : AppleClang v17.0.0.17000013 (/usr/bin/cc) - Fortran : GNU v15.1.0 (/opt/homebrew/bin/gfortran) + C : NVHPC v23.11.0 (/opt/nvidia/hpc_sdk/Linux_x86_64/23.11/compilers/bin/nvc) + Fortran : NVHPC v23.11.0 (/opt/nvidia/hpc_sdk/Linux_x86_64/23.11/compilers/bin/nvfortran) PRE_PROCESS : OFF SIMULATION : OFF @@ -94,31 +94,31 @@ syscheck: OpenACC : OFF OpenMP : OFF - Fypp : /Users/dan/Documents/repos/MFC/build/venv/bin/fypp + Fypp : /home/dan/Documents/repos/MFC/build/venv/bin/fypp Doxygen : Build Type : Release Configuration Environment: - CC : /usr/bin/cc - CXX : /usr/bin/c++ - FC : /opt/homebrew/bin/gfortran + CC : nvc + CXX : nvc++ + FC : nvfortran OMPI_CC : OMPI_CXX : OMPI_FC : -simulation: +pre_process: CMake Configuration: - CMake v4.0.3 on oppenheimer + CMake v3.28.3 on schwarzschild - C : AppleClang v17.0.0.17000013 (/usr/bin/cc) - Fortran : GNU v15.1.0 (/opt/homebrew/bin/gfortran) + C : NVHPC v23.11.0 (/opt/nvidia/hpc_sdk/Linux_x86_64/23.11/compilers/bin/nvc) + Fortran : NVHPC v23.11.0 (/opt/nvidia/hpc_sdk/Linux_x86_64/23.11/compilers/bin/nvfortran) - PRE_PROCESS : OFF - SIMULATION : ON + PRE_PROCESS : ON + SIMULATION : OFF POST_PROCESS : OFF SYSCHECK : OFF DOCUMENTATION : OFF @@ -128,16 +128,16 @@ simulation: OpenACC : OFF OpenMP : OFF - Fypp : /Users/dan/Documents/repos/MFC/build/venv/bin/fypp + Fypp : /home/dan/Documents/repos/MFC/build/venv/bin/fypp Doxygen : Build Type : Release Configuration Environment: - CC : /usr/bin/cc - CXX : /usr/bin/c++ - FC : /opt/homebrew/bin/gfortran + CC : nvc + CXX : nvc++ + FC : nvfortran OMPI_CC : OMPI_CXX : OMPI_FC : @@ -145,10 +145,50 @@ simulation: CPU: CPU Info: - From sysctl -a - machdep.cpu.cores_per_package: 8 - machdep.cpu.core_count: 8 - machdep.cpu.logical_per_package: 8 - machdep.cpu.thread_count: 8 - machdep.cpu.brand_string: Apple M2 + From lscpu + Architecture: x86_64 + CPU op-mode(s): 32-bit, 64-bit + Address sizes: 46 bits physical, 48 bits virtual + Byte Order: Little Endian + CPU(s): 20 + On-line CPU(s) list: 0-19 + Vendor ID: GenuineIntel + Model name: 12th Gen Intel(R) Core(TM) i7-12700K + CPU family: 6 + Model: 151 + Thread(s) per core: 2 + Core(s) per socket: 12 + Socket(s): 1 + Stepping: 2 + CPU(s) scaling MHz: 34% + CPU max MHz: 5100.0000 + CPU min MHz: 800.0000 + BogoMIPS: 7219.20 + Flags: fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush dts acpi mmx fxsr sse sse2 ss ht tm pbe syscall nx pdpe1gb rdtscp lm constant_tsc art arch_perfmon pebs bts rep_good nopl xtopology nonstop_tsc cpuid aperfmperf tsc_known_freq pni pclmulqdq dtes64 monitor ds_cpl vmx smx est tm2 ssse3 sdbg fma cx16 xtpr pdcm pcid sse4_1 sse4_2 x2apic movbe popcnt tsc_deadline_timer aes xsave avx f16c rdrand lahf_lm abm 3dnowprefetch cpuid_fault cat_l2 cdp_l2 ssbd ibrs ibpb stibp ibrs_enhanced tpr_shadow flexpriority ept vpid ept_ad fsgsbase tsc_adjust bmi1 avx2 smep bmi2 erms invpcid rdt_a rdseed adx smap clflushopt clwb intel_pt sha_ni xsaveopt xsavec xgetbv1 xsaves split_lock_detect user_shstk avx_vnni dtherm ida arat pln pts hwp hwp_notify hwp_act_window hwp_epp hwp_pkg_req hfi vnmi umip pku ospke waitpkg gfni vaes vpclmulqdq rdpid movdiri movdir64b fsrm md_clear serialize pconfig arch_lbr ibt flush_l1d arch_capabilities + Virtualization: VT-x + L1d cache: 512 KiB (12 instances) + L1i cache: 512 KiB (12 instances) + L2 cache: 12 MiB (9 instances) + L3 cache: 25 MiB (1 instance) + NUMA node(s): 1 + NUMA node0 CPU(s): 0-19 + Vulnerability Gather data sampling: Not affected + Vulnerability Ghostwrite: Not affected + Vulnerability Indirect target selection: 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 Old microcode: Not affected + Vulnerability Reg file data sampling: Mitigation; Clear Register File + Vulnerability Retbleed: Not affected + Vulnerability Spec rstack overflow: Not affected + 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; Enhanced / Automatic IBRS; IBPB conditional; PBRSB-eIBRS SW sequence; BHI BHI_DIS_S + Vulnerability Srbds: Not affected + Vulnerability Tsa: Not affected + Vulnerability Tsx async abort: Not affected + Vulnerability Vmscape: Mitigation; IBPB before exit to userspace diff --git a/tests/127A967A/golden.txt b/tests/127A967A/golden.txt index 461c051fce..c03db9cc83 100644 --- a/tests/127A967A/golden.txt +++ b/tests/127A967A/golden.txt @@ -1,10 +1,10 @@ D/cons.1.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.1.00.000050.dat 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000005 1.00000000000043 1.0000000000026 1.00000000005959 1.0000000024797 1.00000002884774 1.000000268344 1.0000022492814 1.00000881795754 1.00001371840713 1.00000443840522 1.00000090815759 1.00000002163506 0.99999999546873 0.99999999893431 0.99999999995108 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000003 1.00000000000025 1.00000000000164 1.00000000006741 1.00000000310871 1.00000003909965 1.0000004222109 1.00000388996087 1.00002981279292 1.0001121495435 1.00015970597748 1.00005030821396 1.00000704532692 0.9999997784402 0.99999986746852 0.99999997887725 0.99999999810318 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000001 1.00000000000013 1.00000000000067 1.00000000004847 1.00000000280479 1.00000003912007 1.00000047005507 1.00000481915722 1.00004304365093 1.00028137348933 1.00098398898122 1.00119247216308 1.00033155648737 1.00002815948502 0.9999925925958 0.99999778158853 0.99999965418174 0.99999997292247 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000001 1.00000000000006 1.00000000000032 1.00000000002132 1.00000000196692 1.00000003143912 1.00000041066882 1.00000464290624 1.00004509739957 1.00039546155216 1.00236592891775 1.00832050844317 1.0079562335406 1.00121590984065 0.99990565023817 0.9998904135531 0.999973145391 0.99999579630136 0.99999968060331 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000006 1.00000000000036 1.00000000001091 1.00000000128027 1.00000002163753 1.00000030956655 1.00000388049165 1.00004125758829 1.00034805294784 1.00309002680186 1.01219267889909 1.02200070948297 1.02056773477218 1.00153837617589 0.99651380512362 0.99865459424484 0.99972277713394 0.9999679841293 0.99999724213128 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000001 1.00000000000006 1.00000000000021 1.00000000037969 1.00000000754417 1.00000011387846 1.00000151935601 1.00001734990726 1.0001660790427 1.00126282138318 1.0108139412836 1.03206067147984 1.05300862232409 1.01339376889326 0.98106287485993 0.98433200323471 0.99433836821428 0.99898659544772 0.99986984774938 0.99998812397824 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000001 1.00000000000008 1.00000000000506 1.00000000123552 1.00000002035061 1.00000029939552 1.00000381817864 1.00004112973825 1.0003659128127 1.00269997063055 1.01799970198972 1.04080379540375 1.04076142972022 0.97385083642864 0.97350359935548 0.96732167119668 0.98538130677003 0.99673543657761 0.99961347529604 0.99996482794384 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000001 1.00000000000005 1.0000000000002 1.00000000040986 1.00000000791219 1.00000011686761 1.00000151868187 1.00001681362642 1.00015455973484 1.00121854100666 1.00740116753425 1.02043122657822 1.02017143651674 0.97720805419775 0.95062397167624 0.97342050771063 0.99048814308727 0.99867471003173 0.99982010941489 0.99998402836755 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000006 1.00000000000039 1.00000000001762 1.00000000153182 1.00000002462995 1.00000034103265 1.00000407572377 1.00004028342167 1.00033737345661 1.00179209278649 1.00523046934847 0.99863878548162 0.98057005063137 0.98029529880362 0.9892094248511 0.9968940645618 0.999617805115 0.99995280635024 0.99999606288215 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000001 1.00000000000008 1.0000000000004 1.00000000003114 1.00000000221573 1.00000003303275 1.00000040084223 1.00000415882374 1.00003670896603 1.00017924492199 1.00034182020691 0.99889608004209 0.9913663848813 0.99093596835175 0.99737075453252 0.99953102233892 0.99994459312503 0.99999413957337 0.99999953734667 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000002 1.00000000000016 1.00000000000068 1.00000000007248 1.00000000334943 1.00000004396931 1.00000051603952 1.00000366201079 1.00001426891511 0.99999294749085 0.99967438316695 0.99839194696475 0.99875783006288 0.99966333748857 0.99994558580123 0.99999361788977 0.99999935684644 0.99999995079931 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000004 1.0000000000003 1.00000000000245 1.00000000007772 1.00000000330446 1.00000003885607 1.000000270074 1.0000008055191 0.99999429022657 0.99993193529656 0.99978764529882 0.99985323288181 0.99996081194377 0.99999460413534 0.99999939079677 0.99999994179094 0.99999999571316 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000001 1.00000000000007 1.00000000000057 1.00000000000426 1.00000000011979 1.00000000233761 1.00000001579327 1.00000002119145 0.99999910269408 0.99999279365215 0.99997513392116 0.99998497486989 0.99999626159834 0.999999542977 0.99999994999873 0.99999999539704 0.99999999984392 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000001 1.00000000000007 1.00000000000059 1.00000000000434 1.00000000007903 1.00000000075092 0.99999999785989 0.99999988352251 0.99999929692791 0.99999760923474 0.99999869190155 0.99999969642163 0.99999996585307 0.99999999639562 0.99999999985522 0.99999999999626 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000005 1.00000000000037 1.00000000000063 0.99999999999913 0.99999999970334 0.99999999081601 0.99999994209594 0.99999980480045 0.99999990110823 0.99999997848403 0.9999999978239 0.99999999991293 0.999999999998 0.99999999999941 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000004 1.00000000000004 0.99999999999921 0.99999999998809 0.99999999947736 0.99999999578767 0.9999999860878 0.99999999333684 0.99999999868546 0.99999999996919 0.99999999999952 0.9999999999997 0.99999999999992 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.99999999999991 0.99999999999914 0.99999999999596 0.99999999981953 0.99999999915024 0.99999999967212 0.9999999999897 0.9999999999998 0.99999999999988 0.99999999999995 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.99999999999989 0.99999999999949 0.99999999999962 0.99999999999994 0.99999999999987 0.99999999999968 0.99999999999996 0.99999999999998 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.99999999999999 0.9999999999999 0.9999999999999 0.99999999999997 0.99999999999996 0.99999999999996 0.99999999999999 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.99999999999999 0.99999999999998 0.99999999999999 0.99999999999999 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 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.1.00.000050.dat 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000006 1.00000000000043 1.00000000000288 1.00000000007181 1.00000000271538 1.00000003114292 1.00000028271079 1.00000237810608 1.00000947595553 1.00001669863406 1.00000578523373 1.00000128163561 1.00000005052877 0.99999999818657 0.99999999918809 0.9999999999607 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000004 1.00000000000026 1.00000000000175 1.00000000008057 1.00000000336931 1.00000004194509 1.00000044804715 1.00000398672278 1.00003113000359 1.00011849518368 1.00019281640024 1.00007539339182 1.00001029284701 1.00000008513262 0.99999990641924 0.99999998026692 0.99999999803018 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000001 1.00000000000013 1.00000000000065 1.000000000059 1.00000000303853 1.00000004184774 1.00000049758057 1.00000500473337 1.00004295498666 1.00029051852944 1.00102617252298 1.0015578066129 1.00039245399241 1.0000446849798 0.9999950492734 0.99999816279192 0.99999970895921 0.99999997655521 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000001 1.00000000000006 1.00000000000031 1.00000000002764 1.00000000214045 1.00000003370695 1.00000043728209 1.00000488670849 1.0000466042131 1.00039960917671 1.00237890591773 1.00817369176186 1.00742475396269 1.00112891146726 0.99995231743186 0.99990419348023 0.99997551948358 0.9999961668255 0.99999970382589 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000006 1.00000000000037 1.0000000000137 1.00000000138932 1.00000002312672 1.0000003279431 1.00000406513344 1.00004265400048 1.00035428231758 1.00306535127726 1.01176432964839 1.02080849461241 1.01849363166686 1.00071967851098 0.99612715569675 0.99867228921986 0.99972566039978 0.99996901027178 0.99999734974589 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000001 1.00000000000006 1.00000000000027 1.00000000041769 1.000000008111 1.00000012137906 1.00000160108769 1.0000180291346 1.00016982000165 1.00126776460281 1.01056367892985 1.03111573459385 1.05045441740556 1.0095179590934 0.97963160144655 0.98437220106963 0.9945921123739 0.99907244030211 0.99987858359002 0.99998893583312 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000001 1.00000000000006 1.00000000000972 1.00000000137491 1.00000002239943 1.00000032727745 1.00000413818502 1.00004414243358 1.00038742791507 1.00282599396604 1.01735007480092 1.03983534035156 1.03103462270783 0.9968723253116 0.9760539968657 0.96630514763966 0.98546491601524 0.99697061977956 0.99963311052297 0.99996700400834 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000001 1.00000000000007 1.0000000000003 1.0000000004559 1.00000000858716 1.00000012580247 1.00000161782602 1.00001768975576 1.0001601233277 1.0012442046411 1.00722839855131 1.01991493114609 1.01854438892549 0.98064565129951 0.95266468820734 0.97270369211447 0.99049037520195 0.99874172083673 0.99982835241773 0.99998480985491 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000006 1.00000000000036 1.00000000002295 1.0000000016717 1.00000002660311 1.00000036672983 1.00000435564621 1.00004259758901 1.00035042259261 1.00181367237044 1.00496161581888 0.99902937064607 0.98012910226079 0.98075426412036 0.98913020423382 0.99695101388679 0.99963546135761 0.99995492655368 0.99999624684937 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000001 1.00000000000008 1.00000000000034 1.00000000004051 1.00000000244468 1.00000003597809 1.00000043537477 1.00000451163027 1.00003902194331 1.00018466296273 1.00035208350989 0.99887347256755 0.991531199359 0.99107828005894 0.99741530988574 0.99954147598881 0.99994638794922 0.99999437190266 0.99999955838756 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000002 1.00000000000016 1.0000000000006 1.00000000009411 1.00000000375863 1.00000004915902 1.0000005784603 1.00000404252199 1.00001579977966 0.99999914048408 0.99967311517969 0.99861215177946 0.99878600061569 0.99966967295942 0.99994636257586 0.99999380628784 0.99999938120808 0.99999995290285 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000006 1.00000000000038 1.00000000000255 1.00000000011006 1.00000000390659 1.00000004604577 1.00000031723401 1.00000103075958 0.99999550941245 0.99993300708815 0.99980319434907 0.99985726584129 0.99996164866634 0.99999470913309 0.99999941081333 0.99999994405404 0.99999999589962 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000007 1.00000000000045 1.0000000000027 1.00000000009833 1.00000000296202 1.00000001980138 1.00000003930142 0.99999910436678 0.99999312942148 0.99997696655411 0.99998543996395 0.99999634798387 0.99999955484211 0.99999995183479 0.99999999558412 0.99999999985705 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000001 1.00000000000008 1.00000000000057 1.00000000000503 1.00000000011666 1.00000000102043 0.99999999927868 0.99999988340924 0.99999934236158 0.99999777188146 0.99999873643365 0.99999970433743 0.99999996698092 0.99999999654613 0.99999999986614 0.99999999999645 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000006 1.00000000000039 1.00000000000098 1.00000000000198 0.99999999984365 0.99999999096289 0.99999994640167 0.99999981771467 0.99999990478351 0.99999997912095 0.99999999791682 0.9999999999206 0.99999999999806 0.99999999999943 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000004 1.00000000000017 0.9999999999999 0.99999999999938 0.99999999949562 0.99999999613266 0.99999998699973 0.99999999360585 0.99999999873218 0.99999999997249 0.99999999999952 0.9999999999997 0.99999999999992 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000002 0.99999999999995 0.99999999999973 0.99999999999668 0.99999999984281 0.99999999921433 0.99999999969006 0.99999999999088 0.99999999999979 0.99999999999988 0.99999999999995 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.99999999999995 0.99999999999955 0.99999999999958 0.99999999999994 0.99999999999987 0.99999999999971 0.99999999999996 0.99999999999998 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.99999999999999 0.99999999999991 0.99999999999989 0.99999999999996 0.99999999999996 0.99999999999996 0.99999999999999 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.99999999999999 0.99999999999998 0.99999999999999 0.99999999999999 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 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.2.00.000000.dat 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 -D/cons.2.00.000050.dat 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.04999999999998 0.04999999999985 0.04999999999852 0.04999999995664 0.0499999984001 0.04999997957838 0.04999978232802 0.04999797157363 0.0499902878356 0.04998182739876 0.0499960751942 0.04999938106452 0.05000003381 0.05000000856874 0.05000000100677 0.05000000001495 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.04999999999998 0.04999999999991 0.04999999999924 0.04999999996505 0.04999999804289 0.0499999723654 0.04999967362416 0.04999666110062 0.04997179167587 0.04987362208779 0.04979500270019 0.04995332370807 0.04999590208714 0.05000077234015 0.0500001640228 0.05000002054782 0.05000000142656 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.04999999999997 0.04999999999991 0.04999999998707 0.0499999984459 0.04999997630444 0.04999968978865 0.04999651432822 0.04996513200567 0.04973947300454 0.04886228959382 0.04844702739997 0.04965988470099 0.05002665801026 0.05001300315021 0.05000224611704 0.05000027001187 0.05000001833414 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.04999999999997 0.04999999999993 0.04999999999829 0.04999999908893 0.04999998433892 0.04999977655061 0.04999724417788 0.04997068402071 0.04971003845997 0.04791717459303 0.04046153933614 0.03920886810765 0.04958772313699 0.05050682767823 0.0501543137144 0.05002283678706 0.05000246455708 0.05000017634278 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.04999999999999 0.04999999999998 0.04999999999171 0.04999999947224 0.04999999121361 0.04999986161475 0.04999804984207 0.04997631063731 0.04976510179399 0.04756388821286 0.03707798526772 0.0151910726345 0.02137024692702 0.05415951024026 0.05639912039192 0.05154645384744 0.05020172938319 0.05001953443281 0.05000145015627 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.04999999999994 0.04999999987455 0.04999999765034 0.04999996265127 0.04999945547732 0.04999314818286 0.04992741273364 0.04937334916842 0.04445347192799 0.0277792189321 0.01595995610544 0.0044145149174 0.04654566530577 0.06866919011405 0.0553730268946 0.05072330212586 0.05007789766397 0.05000600781569 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05000000000024 0.05000000005581 0.05000000091144 0.05000001339099 0.05000017175019 0.05000187595906 0.05001720123517 0.05012957278729 0.05111144878045 0.05086645013896 0.00453373308381 0.00424225920536 0.00424074658185 0.05235222860661 0.05013473202561 0.04995480054216 0.04999204966421 0.04999904415738 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05000000000008 0.05000000016857 0.05000000314246 0.05000004900437 0.05000069772611 0.05000858056632 0.05008923955642 0.05075540985388 0.05684935329741 0.07583370956532 0.05071451110407 0.00425688381464 0.02115766088615 0.032375754643 0.04423933031983 0.04912086777526 0.04990489321638 0.04999248838338 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05000000000001 0.05000000000005 0.05000000001203 0.05000000074876 0.0500000123787 0.05000018606488 0.05000248115519 0.05002811163529 0.05027191286947 0.05210260596651 0.05973099420368 0.05363544462246 0.02150507607173 0.0179680479152 0.03805057168408 0.04724005011161 0.04970675862522 0.04996812419348 0.04999761281952 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05000000000001 0.05000000000005 0.05000000000011 0.05000000001001 0.05000000129132 0.05000002119048 0.05000028703293 0.05000331511974 0.05003487805975 0.05023205427589 0.05094292163251 0.04992463698527 0.03682876661385 0.03906689587369 0.04748566506318 0.04960820666813 0.04995797035164 0.04999596878183 0.0499997092499 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05000000000001 0.05000000000006 0.05000000000015 0.05000000002887 0.05000000219753 0.05000003129742 0.05000038395344 0.0500037547041 0.05002254076395 0.05006241989968 0.04969127516492 0.04771073658908 0.0484425882517 0.049663299753 0.04995030009267 0.04999471429796 0.04999951342255 0.04999996581089 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05000000000003 0.05000000000013 0.05000000000127 0.05000000005615 0.05000000265182 0.05000003593323 0.05000031670824 0.05000161582414 0.05000070415837 0.0499303255639 0.04970371214049 0.04982060256319 0.04995961525821 0.04999483410635 0.049999469564 0.04999995336772 0.04999999687562 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05000000000004 0.05000000000028 0.05000000000267 0.05000000008517 0.05000000248336 0.05000002113343 0.05000008511463 0.04999943086575 0.04999255070326 0.04996587640856 0.04998184990993 0.04999606952128 0.04999954815481 0.04999995490458 0.04999999626085 0.04999999988266 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05000000000006 0.05000000000056 0.05000000000223 0.050000000062 0.05000000112714 0.05000000268971 0.04999989964493 0.04999925445991 0.04999676400779 0.04999843069989 0.04999967611215 0.04999996550212 0.04999999673104 0.04999999985583 0.04999999999802 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05000000000006 0.05000000000045 0.0500000000012 0.04999999999997 0.0499999997887 0.04999999173111 0.04999993696475 0.04999973869117 0.04999988195783 0.0499999768544 0.04999999774288 0.0499999998974 0.04999999999869 0.04999999999946 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05000000000004 0.05000000000017 0.04999999999968 0.04999999998617 0.04999999948497 0.04999999540003 0.04999998160065 0.04999999212809 0.04999999855683 0.04999999996031 0.04999999999938 0.04999999999966 0.04999999999992 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05000000000001 0.04999999999995 0.04999999999907 0.04999999999818 0.04999999980571 0.04999999889115 0.04999999960881 0.04999999999303 0.04999999999979 0.04999999999987 0.04999999999995 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.04999999999988 0.04999999999938 0.04999999999949 0.04999999999973 0.04999999999979 0.04999999999962 0.04999999999996 0.04999999999998 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.04999999999999 0.04999999999989 0.04999999999988 0.04999999999996 0.04999999999995 0.04999999999995 0.04999999999999 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.04999999999999 0.04999999999998 0.04999999999999 0.04999999999999 0.04999999999999 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 +D/cons.2.00.000050.dat 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.04999999999998 0.04999999999984 0.0499999999983 0.04999999994814 0.04999999827223 0.04999997809365 0.04999977103331 0.04999789084471 0.04998980785733 0.0499777431431 0.04999478207362 0.04999901717531 0.05000001468687 0.05000000688529 0.0500000009307 0.05000000001235 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.04999999999998 0.04999999999991 0.04999999999913 0.04999999995849 0.04999999788855 0.04999997039532 0.04999965329071 0.04999654967292 0.04997096476554 0.04986926191391 0.04975157073011 0.04992535781146 0.04999299168076 0.05000056321692 0.05000013956555 0.05000001818762 0.05000000131896 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.04999999999997 0.0499999999999 0.04999999998288 0.04999999832456 0.04999997459636 0.04999966986964 0.04999633071829 0.04996441239606 0.04973713064698 0.04885214899022 0.04795430514611 0.04961178123341 0.05002104702808 0.05001124074526 0.0500020493349 0.05000024615412 0.05000001683448 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.04999999999996 0.04999999999993 0.04999999999705 0.04999999900718 0.049999983113 0.04999976051955 0.04999707024754 0.04996926393913 0.04970473592139 0.04795126004587 0.04073248544198 0.03870144311069 0.04998313108272 0.05050556662046 0.05014421575785 0.05002192150746 0.05000233218748 0.05000016719557 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.04999999999999 0.04999999999997 0.04999999999054 0.04999999942748 0.04999999056475 0.04999985244534 0.04999794613869 0.04997541739052 0.04976070579744 0.04754837644817 0.03750705642905 0.01687610131467 0.0260104744552 0.05483656292807 0.05648414360292 0.05146983521564 0.05019013600146 0.05001856395631 0.05000137037869 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.04999999999986 0.04999999985719 0.04999999739493 0.04999995877253 0.04999940324412 0.04999254946397 0.04992178068844 0.04932730921101 0.04450861311279 0.0284064892654 0.01810998989882 0.00408933428452 0.04582939827219 0.06900731097693 0.05511108475807 0.05065740619744 0.05007254216685 0.05000558701133 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05000000000046 0.05000000006161 0.05000000099502 0.05000001448437 0.05000018391705 0.05000198513085 0.05001790544163 0.05013210529804 0.05106399001095 0.050778977922 0.00417649353653 0.00403810961506 0.00395377916458 0.05251834756717 0.05019101501505 0.04996889604483 0.04999306351946 0.04999914140432 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05000000000017 0.05000000019067 0.050000003464 0.05000005378034 0.05000075937697 0.05000924895569 0.05009509971706 0.05079891156805 0.05675360711629 0.07516163598912 0.04774417610564 0.00397237894255 0.02178505731377 0.03182422019579 0.04445846709227 0.04918924904851 0.04991066893273 0.04999297791222 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05000000000001 0.05000000000007 0.05000000001424 0.05000000081949 0.05000001340004 0.05000020045914 0.05000265223812 0.05002973267657 0.05028372717268 0.05214717965797 0.05939858565882 0.05443889480257 0.01925948160857 0.01866860332088 0.03789932884247 0.04732369276977 0.04972147974611 0.04996963089153 0.04999772628695 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05000000000001 0.05000000000005 0.0500000000001 0.05000000001325 0.05000000141565 0.05000002303351 0.05000031031362 0.0500035687506 0.05003721678124 0.05024147035325 0.05092153964575 0.04983342176267 0.03773030319965 0.03920739752237 0.04751280031407 0.04962141101633 0.04995974085185 0.04999615651886 0.04999972421136 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05000000000001 0.05000000000005 0.05000000000014 0.05000000003803 0.05000000242311 0.0500000343417 0.05000042167454 0.05000407508248 0.05002408997597 0.05006363430152 0.04964286829612 0.04804151142116 0.04846301476025 0.04966702968453 0.04995168179628 0.0499949150914 0.04999953441045 0.04999996740659 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05000000000003 0.05000000000014 0.05000000000152 0.05000000007407 0.05000000303725 0.05000004079769 0.05000035598793 0.05000183033478 0.05000205220715 0.04993010966409 0.04972621373822 0.04982439766471 0.0499602599335 0.04999497100147 0.04999948904611 0.0499999552915 0.04999999701598 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05000000000005 0.05000000000036 0.050000000003 0.05000000009808 0.05000000297638 0.05000002495788 0.05000010699889 0.04999936842213 0.04999287570584 0.04996842732521 0.04998231653813 0.04999614407104 0.0499995616287 0.04999995661285 0.04999999641878 0.0499999998931 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05000000000005 0.05000000000047 0.05000000000164 0.0500000000961 0.05000000142397 0.05000000415143 0.04999989510804 0.04999930561739 0.04999698548698 0.04999847748197 0.04999968335433 0.04999996669589 0.04999999686299 0.04999999986651 0.04999999999806 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05000000000008 0.05000000000045 0.05000000000144 0.05000000000298 0.04999999999738 0.04999999159285 0.04999994199783 0.04999975603595 0.04999988593342 0.04999997746177 0.04999999783465 0.04999999990586 0.04999999999864 0.04999999999947 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05000000000004 0.05000000000027 0.05 0.04999999999931 0.04999999948895 0.04999999580634 0.04999998281984 0.0499999924225 0.04999999860067 0.04999999996471 0.04999999999939 0.04999999999966 0.04999999999993 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05000000000003 0.04999999999998 0.04999999999987 0.04999999999839 0.04999999983369 0.04999999897414 0.04999999962907 0.04999999999417 0.04999999999977 0.04999999999987 0.04999999999995 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.04999999999995 0.04999999999944 0.04999999999945 0.04999999999974 0.0499999999998 0.04999999999965 0.04999999999996 0.04999999999998 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.0499999999999 0.04999999999987 0.04999999999996 0.04999999999995 0.04999999999995 0.04999999999999 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.04999999999999 0.04999999999998 0.04999999999999 0.04999999999999 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 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 -5e-14 -3.9e-13 -1.67e-12 -3.223e-11 -1.59892e-09 -1.681843e-08 -1.3856516e-07 -1.19948516e-06 -2.46667983e-06 4.7305218e-07 2.66766224e-06 6.2837466e-07 5.273083e-08 1.96152e-09 -5.4952e-10 -4.108e-11 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -3e-14 -2.6e-13 -1.37e-12 -5.432e-11 -2.16736e-09 -2.660954e-08 -2.6697658e-07 -2.18536472e-06 -1.746985815e-05 -3.188577043e-05 7.52103543e-06 3.632087106e-05 7.51955303e-06 4.9759222e-07 -1.338751e-08 -1.069859e-08 -9.4547e-10 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -2e-14 -1.4e-13 -7.6e-13 -5.059e-11 -2.21996e-09 -3.056143e-08 -3.5028609e-07 -3.39228723e-06 -2.746387943e-05 -0.00018575429547 -0.000274934958 9.539710802e-05 0.00033005014801 6.675412011e-05 3.30020203e-06 -1.11297457e-06 -2.1984161e-07 -1.579524e-08 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -1e-14 -6e-14 -3.4e-13 -2.554e-11 -1.82152e-09 -2.788533e-08 -3.5489957e-07 -3.90037565e-06 -3.648087169e-05 -0.00029578458508 -0.00179743546873 -0.00239494487606 0.00161491683665 0.0024996647428 0.00038130553297 5.24133202e-06 -2.007014129e-05 -3.61077303e-06 -2.4754039e-07 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -1e-14 -6e-14 -4.5e-13 -6.31e-12 -1.30044e-09 -2.123641e-08 -3.0063779e-07 -3.68262057e-06 -3.795807272e-05 -0.00030550485877 -0.00263302306327 -0.0090384013136 -0.00873250320269 0.0078136994975 0.0124105086927 0.00114908348924 -0.00068760738192 -0.00022594914669 -2.784193619e-05 -2.30527625e-06 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -1e-14 -7e-14 -3e-13 -4.1959e-10 -8.25765e-09 -1.2530581e-07 -1.66468922e-06 -1.893395559e-05 -0.0001808096265 -0.00137359789804 -0.01204430745508 -0.03908406264867 -0.05057186616935 -0.09237481699761 -0.0321005396525 -0.01189143877782 -0.00562583379135 -0.00093424222551 -0.00012897177645 -1.135233312e-05 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -1e-14 -9e-14 -6.33e-12 -1.53255e-09 -2.588004e-08 -3.8428663e-07 -4.95617336e-06 -5.421356328e-05 -0.00049271119826 -0.00374025013747 -0.02578374947252 -0.06762679913976 -0.09486948662964 -0.0887703235992 -0.08873867157797 -0.04284919020004 -0.02294038586894 -0.00448841560301 -0.00052374594717 -4.868601861e-05 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -1e-14 -7e-14 -2.8e-13 -4.5901e-10 -8.75455e-09 -1.295397e-07 -1.67007906e-06 -1.831677152e-05 -0.00016661899785 -0.00129817882927 -0.00820633189809 -0.01762033243704 -0.03822939550003 -0.0890763471673 -0.04796540451312 -0.03315006370455 -0.01044006199519 -0.00142376734459 -0.00019910533328 -1.732140606e-05 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -1e-14 -6e-14 -4.8e-13 -1.282e-11 -1.56065e-09 -2.477109e-08 -3.3804155e-07 -3.92518345e-06 -3.702744707e-05 -0.000292721984 -0.00114497263143 0.0009781418049 0.01226414392766 0.00533697016287 -0.00887069847331 -0.00866653199124 -0.00265194037207 -0.00033177136859 -4.38380998e-05 -3.48713425e-06 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -1e-14 -8e-14 -4.6e-13 -3.362e-11 -2.00578e-09 -2.835812e-08 -3.2349586e-07 -3.07803939e-06 -2.633475127e-05 -5.532420641e-05 0.00042317770559 0.00309383506817 0.00194882636492 -0.00281212147835 -0.00209095279699 -0.00035389951338 -4.519646062e-05 -5.05469823e-06 -3.7183709e-07 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -3e-14 -1.7e-13 -8.1e-13 -7.518e-11 -2.69171e-09 -3.522605e-08 -4.029288e-07 -2.31994601e-06 2.15440461e-06 7.707780667e-05 0.0004411811563 0.00013708614344 -0.00037533011449 -0.00023703056799 -3.54902698e-05 -4.56482649e-06 -4.9357931e-07 -3.451015e-08 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -4e-14 -2.9e-13 -2.1e-12 -5.694e-11 -2.08423e-09 -2.255534e-08 -1.1020425e-07 5.1510156e-07 9.69046654e-06 5.225840647e-05 1.047315795e-05 -4.506315655e-05 -2.41701399e-05 -3.11726451e-06 -3.949096e-07 -4.089622e-08 -2.67373e-09 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -1e-14 -7e-14 -5.2e-13 -3.25e-12 -9.159e-11 -1.36576e-09 -1.13199e-09 5.760326e-08 9.2156428e-07 4.98475148e-06 9.0767276e-07 -4.53853721e-06 -2.0576803e-06 -2.3981737e-07 -2.989892e-08 -2.91138e-09 -9.446e-11 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -1e-14 -6e-14 -3e-13 -3.24e-12 -6.4e-11 1.0862e-10 5.9088e-09 7.972211e-08 4.0379866e-07 6.45278e-08 -3.8489106e-07 -1.5048811e-07 -1.652856e-08 -2.0011e-09 -6.705e-11 -3.05e-12 -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 -1e-13 1.4e-13 1.59e-12 2.8808e-10 5.51805e-09 2.949187e-08 3.82064e-09 -2.811843e-08 -9.95829e-09 -1.0167e-09 -2.363e-11 -1.55e-12 -2.6e-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 -1e-14 1.1e-13 9.2e-13 1.43e-12 2.8863e-10 2.03259e-09 1.7895e-10 -1.91178e-09 -5.8964e-10 -3.85e-12 -7e-14 -9e-14 -2e-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 1e-13 2.9e-13 3.97e-12 8.724e-11 1.343e-11 -9.632e-11 -8.05e-12 -6e-14 -2e-14 -2e-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 7e-14 1e-13 1.7e-13 4e-14 -2.9e-13 -7e-14 -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 1e-14 1e-14 -1e-14 -0.0 -0.0 -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 +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 -5e-14 -3.9e-13 -1.78e-12 -3.911e-11 -1.74922e-09 -1.833428e-08 -1.4739532e-07 -1.27736878e-06 -3.04747498e-06 4.6766061e-07 3.15825173e-06 7.9190195e-07 7.130525e-08 3.39622e-09 -1.9677e-10 -3.282e-11 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -3e-14 -2.6e-13 -1.39e-12 -6.535e-11 -2.34913e-09 -2.865759e-08 -2.8435608e-07 -2.23793329e-06 -1.844429886e-05 -3.97166403e-05 6.23584408e-06 4.396759076e-05 9.79356726e-06 7.09607e-07 2.001036e-08 -1.134446e-08 -1.2181e-09 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -2e-14 -1.4e-13 -7.3e-13 -6.065e-11 -2.39532e-09 -3.269023e-08 -3.7002402e-07 -3.49137373e-06 -2.672327436e-05 -0.00019566105515 -0.00036860179297 8.982642397e-05 0.00042025278658 8.289255333e-05 4.75501042e-06 -6.6392712e-07 -1.6551957e-07 -1.217153e-08 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -1e-14 -5e-14 -3.4e-13 -3.243e-11 -1.96465e-09 -2.985791e-08 -3.7683021e-07 -4.08562589e-06 -3.74368336e-05 -0.0002972800528 -0.001852946704 -0.00256540068359 0.00160436005203 0.00264890200916 0.00045576583904 1.437024017e-05 -1.710990521e-05 -3.22548369e-06 -2.2654359e-07 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -1e-14 -6e-14 -4.6e-13 -8.72e-12 -1.39949e-09 -2.265254e-08 -3.1748646e-07 -3.84222098e-06 -3.901391428e-05 -0.00030777687094 -0.00259976268164 -0.00862567627009 -0.00772046069777 0.0085009404957 0.01080348621802 0.00074577204023 -0.0007821638069 -0.00022947630102 -2.745549866e-05 -2.25314612e-06 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -1e-14 -7e-14 -3e-13 -4.6037e-10 -8.84178e-09 -1.3285217e-07 -1.74300457e-06 -1.952046777e-05 -0.00018309684178 -0.00136061105977 -0.0117416231212 -0.03759407876064 -0.04937292449074 -0.09235492229362 -0.03361156921004 -0.0117265297502 -0.00531630478294 -0.00084874580922 -0.00011968006235 -1.054401454e-05 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -1e-14 -7e-14 -1.203e-11 -1.704e-09 -2.853223e-08 -4.2074747e-07 -5.38259741e-06 -5.834467267e-05 -0.00052374328653 -0.00392894833971 -0.02533875493608 -0.06389830052376 -0.09432335661242 -0.09119804686139 -0.08929349916257 -0.04355499038019 -0.02249449099388 -0.00416113319882 -0.00049638696605 -4.562774963e-05 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -2e-14 -9e-14 -3.6e-13 -5.1019e-10 -9.48255e-09 -1.3901803e-07 -1.7718862e-06 -1.916631173e-05 -0.00017133354094 -0.00131204992887 -0.0080100599856 -0.0171180800931 -0.03808364905326 -0.08971356290152 -0.04636710826939 -0.03366605317458 -0.01058246929373 -0.00136357087711 -0.0001906226981 -1.65495184e-05 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -1e-14 -6e-14 -4.4e-13 -1.801e-11 -1.68987e-09 -2.670771e-08 -3.6290309e-07 -4.19075426e-06 -3.911233508e-05 -0.0003033629295 -0.00115865186389 0.00104298336409 0.0125143712228 0.00502294665182 -0.00887884733542 -0.00859823135586 -0.00264602969882 -0.0003191571288 -4.203604469e-05 -3.33783016e-06 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -1e-14 -8e-14 -4e-13 -4.367e-11 -2.19916e-09 -3.095326e-08 -3.5244559e-07 -3.35931615e-06 -2.820662694e-05 -5.5288126e-05 0.00039145716171 0.00295904118264 0.00193100087806 -0.00267062747442 -0.0020459960142 -0.00034976785253 -4.408693105e-05 -4.87732606e-06 -3.5700228e-07 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -3e-14 -1.7e-13 -7.1e-13 -9.759e-11 -3.05432e-09 -4.003838e-08 -4.6110196e-07 -2.65337274e-06 1.50520799e-06 6.939021171e-05 0.00039394254895 0.00014587261546 -0.00033429811312 -0.00023055876351 -3.560924361e-05 -4.46905295e-06 -4.7709083e-07 -3.318502e-08 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -6e-14 -4e-13 -2.09e-12 -8.711e-11 -2.57411e-09 -2.862493e-08 -1.5170179e-07 4.0641462e-07 8.90633204e-06 4.90351963e-05 1.030584011e-05 -4.139942687e-05 -2.354971821e-05 -3.08528394e-06 -3.8349322e-07 -3.940281e-08 -2.56542e-09 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -5e-14 -3.5e-13 -8e-13 -4.422e-11 -1.76247e-09 -4.11474e-09 5.414014e-08 8.7904568e-07 4.62614954e-06 9.0646319e-07 -4.18824986e-06 -2.00526112e-06 -2.3442709e-07 -2.882936e-08 -2.79481e-09 -8.707e-11 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -1e-14 -6e-14 -3.1e-13 -4.41e-12 -9.177e-11 3.993e-11 5.41937e-09 7.659821e-08 3.7261935e-07 6.639099e-08 -3.5662089e-07 -1.4639863e-07 -1.596041e-08 -1.92106e-09 -6.159e-11 -2.84e-12 -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 -1.8e-13 -1.9e-13 4.1e-13 2.3367e-10 5.29006e-09 2.72264e-08 4.0028e-09 -2.608343e-08 -9.67333e-09 -9.7267e-10 -2.053e-11 -1.42e-12 -2.5e-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 -2e-14 4e-14 2.1e-13 5.1e-13 2.6778e-10 1.8905e-09 1.8886e-10 -1.78011e-09 -5.6697e-10 -3.53e-12 -7e-14 -9e-14 -2e-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 1e-14 5e-14 2.7e-13 3.29e-12 7.804e-11 1.429e-11 -8.784e-11 -7.49e-12 -6e-14 -2e-14 -2e-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 5e-14 1.2e-13 1.8e-13 3e-14 -3.1e-13 -6e-14 -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 1e-14 2e-14 -1e-14 -1e-14 -0.0 -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 D/cons.4.00.000000.dat 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 -D/cons.4.00.000050.dat 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000001 2.50125000000019 2.50125000000149 2.50125000000902 2.50125000020633 2.50125000859586 2.50125009990995 2.50125092798522 2.50125776827054 2.50128036655153 2.50129708958141 2.50126533269354 2.50125314647192 2.50125007738619 2.50124998457466 2.50124999632175 2.5012499998296 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000001 2.50125000000012 2.50125000000088 2.50125000000571 2.50125000023407 2.50125001079274 2.50125013543255 2.50125146090047 2.50126344310427 2.50135289962077 2.50163610354759 2.50179860589056 2.5014236886659 2.50127444509615 2.50124926343659 2.50124954449979 2.50124992711244 2.50124999342688 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000005 2.50125000000047 2.50125000000233 2.50125000016886 2.50125000974817 2.50125013568922 2.50125162909082 2.5012666868019 2.50139886037906 2.50222162635358 2.50463861534345 2.5053487507676 2.50239328912011 2.50134986739888 2.50122473399134 2.50124235065826 2.50124880357281 2.50124990617146 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000003 2.5012500000002 2.50125000000111 2.50125000007452 2.50125000684156 2.5012501092209 2.5012514256506 2.50126610663117 2.50140632378477 2.50261948497485 2.50943626831346 2.53005146739478 2.52870573913042 2.50549044286766 2.50094612786944 2.50087435385203 2.50115718626278 2.50123541558689 2.50124889133073 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000002 2.5012500000002 2.50125000000126 2.50125000003776 2.50125000445301 2.50125007528898 2.50125107618295 2.50126347940734 2.50139317001159 2.50245627919361 2.51195940275518 2.54355698064923 2.57786319463593 2.5730014865444 2.50700771146297 2.48944703638311 2.49662488468056 2.50029031600907 2.50113896366904 2.50124042343811 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000004 2.5012500000002 2.50125000000074 2.50125000132217 2.50125002628312 2.50125039654668 2.50125528862459 2.50131036132202 2.50182751369159 2.50564075098523 2.53904839612585 2.61444433253171 2.69066666706893 2.55195835754466 2.43777436395064 2.44785679958712 2.48177910712744 2.49774287304 2.50079856685968 2.50120874962168 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000003 2.50125000000028 2.5012500000177 2.50125000432556 2.50125007124179 2.50125104816922 2.50126336748043 2.50139400272701 2.50253154304987 2.51072514545815 2.56511359666288 2.64961750394873 2.65031375084002 2.41423737507538 2.41334373735842 2.390418609701 2.45073972457766 2.48985912754992 2.49989774064867 2.50112690010189 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000004 2.50125000000019 2.50125000000069 2.50125000144242 2.50125002783603 2.50125041133388 2.50125534837345 2.50130925656435 2.50179529201511 2.50555473341187 2.52761544929204 2.57507258951665 2.57498490631672 2.42539199948546 2.33205692040653 2.40887352247904 2.46786531789593 2.49657338192713 2.50061593937057 2.50119374462536 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000002 2.5012500000002 2.50125000000136 2.50125000006224 2.50125000540062 2.50125008679268 2.50125120249382 2.50126438404256 2.50139235105698 2.50244421006043 2.50763273607904 2.5201640049212 2.49688264557053 2.43308386834447 2.4318443444926 2.4631187724784 2.49026411723765 2.49989843655853 2.50108329283855 2.5012361056965 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000004 2.50125000000029 2.50125000000142 2.50125000010953 2.50125000782114 2.50125011663489 2.50125141680021 2.50126471651051 2.50138018414076 2.50188889250438 2.50249641619703 2.49739782312315 2.4705981402399 2.46917761136296 2.49193975463941 2.49959012150078 2.50105405245559 2.50122929440961 2.50124836675558 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000008 2.50125000000056 2.50125000000239 2.50125000025802 2.50125001182842 2.50125015540014 2.50125182468991 2.50126300024185 2.50130105183225 2.50122847792327 2.50009543120039 2.49551693731696 2.49683025880101 2.50005560250179 2.50105714178051 2.5012274064266 2.50124772543994 2.50124982615106 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000001 2.50125000000015 2.50125000000106 2.5012500000084 2.50125000027915 2.50125001169484 2.50125013774414 2.50125096075648 2.50125289911089 2.50123005858608 2.50100838217907 2.50049243783323 2.50072761576767 2.50111087628943 2.50123086302127 2.50124784203058 2.50124979401072 2.50124998484546 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000002 2.50125000000026 2.50125000000201 2.50125000001496 2.5012500004273 2.50125000830271 2.50125005631218 2.50125007839928 2.50124683207181 2.50122441445761 2.50116129641358 2.50119652422962 2.50123672379515 2.5012483784002 2.50124982280454 2.50124998370906 2.50124999944361 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000002 2.50125000000026 2.50125000000208 2.50125000001513 2.50125000028291 2.50125000268667 2.5012499926477 2.50124958745734 2.50124750285244 2.50124147354123 2.50124534483427 2.50124892166211 2.50124987880481 2.50124998722627 2.50124999948208 2.50124999998688 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000019 2.50125000000132 2.50125000000211 2.50124999999693 2.5012499989515 2.50124996745384 2.50124979425801 2.50124930398169 2.5012496481016 2.50124992356497 2.5012499922734 2.50124999968645 2.50124999999297 2.50124999999791 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000001 2.50125000000014 2.50125000000015 2.50124999999721 2.50124999995766 2.50124999814137 2.50124998503333 2.50124995040657 2.50124997629505 2.50124999532647 2.50124999988714 2.50124999999831 2.50124999999894 2.50124999999972 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000001 2.50125000000002 2.50124999999968 2.50124999999695 2.5012499999859 2.50124999935426 2.50124999696732 2.50124999882923 2.50124999996242 2.50124999999932 2.50124999999957 2.50124999999983 2.50124999999999 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50124999999999 2.5012499999996 2.50124999999819 2.50124999999873 2.50124999999982 2.50124999999959 2.50124999999891 2.50124999999984 2.50124999999994 2.50124999999999 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50124999999995 2.50124999999966 2.50124999999965 2.50124999999988 2.50124999999986 2.50124999999985 2.50124999999998 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50124999999998 2.50124999999994 2.50124999999998 2.50124999999997 2.50124999999999 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 +D/cons.4.00.000050.dat 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000001 2.5012500000002 2.50125000000151 2.50125000000998 2.50125000024867 2.50125000941406 2.50125010786597 2.50125097768632 2.50125821496118 2.50128264476031 2.50130731290864 2.50126998030678 2.50125443498646 2.50125017752191 2.50124999399955 2.50124999720585 2.50124999986313 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000001 2.50125000000012 2.5012500000009 2.50125000000608 2.50125000027979 2.50125001169772 2.50125014528932 2.50125155027822 2.50126377608012 2.50135746703895 2.50165809029196 2.50191231818588 2.50151006764683 2.50128566198346 2.50125032602125 2.50124967955572 2.50124993185432 2.50124999316768 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000005 2.50125000000047 2.50125000000227 2.50125000020551 2.50125001056112 2.50125014514594 2.50125172439976 2.50126732691038 2.50139851413618 2.50225351303133 2.5047857905196 2.50660556122594 2.50260410158809 2.50140740993358 2.50123324110001 2.50124367455067 2.50124899403553 2.50124991880663 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000002 2.5012500000002 2.50125000000109 2.50125000009657 2.50125000744613 2.50125011709316 2.50125151796248 2.50126695094332 2.5014115250269 2.50263373829059 2.50948323372254 2.52954189387495 2.52681907767422 2.50520704037706 2.50110947728019 2.5009220561005 2.50116544655094 2.50123670533038 2.5012489721221 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000002 2.5012500000002 2.50125000000129 2.50125000004745 2.50125000483233 2.50125008046511 2.50125114001883 2.50126412024148 2.50139801129091 2.50247785979633 2.51187163222722 2.54206001103855 2.57365066616183 2.56578171748537 2.50415025949248 2.48809670051644 2.49668261122518 2.50029981987123 2.50114250522057 2.50124079596293 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000004 2.50125000000021 2.50125000000094 2.50125000145425 2.50125002825391 2.50125042259719 2.50125557197439 2.50131270788844 2.50184032278499 2.50565571373997 2.53816386663735 2.61102657754906 2.68155146976873 2.54036575298295 2.43288708951544 2.44802878919119 2.4826482224058 2.49803953134531 2.50082885799204 2.50121156899632 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000003 2.50125000000022 2.50125000003404 2.50125000481354 2.50125007841391 2.50125114577438 2.50126448771999 2.50140454978715 2.50260690903343 2.51116720382916 2.5625847961624 2.64692547968807 2.62549660773929 2.47065615324219 2.41878329711589 2.38691056635056 2.45102392829795 2.49067837163664 2.49996644443663 2.5011345177431 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000004 2.50125000000024 2.50125000000105 2.50125000160462 2.50125003021336 2.50125044283306 2.50125569833768 2.5013123554235 2.50181505456138 2.50564679340674 2.52699258542755 2.57322183335915 2.56928614251019 2.43388448572468 2.33889719945036 2.4064125196896 2.4678867904422 2.49681084098377 2.50064506077923 2.50119650324408 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000002 2.50125000000021 2.50125000000126 2.50125000008103 2.50125000589384 2.50125009374794 2.50125129312161 2.50126537198153 2.50140052923387 2.50249047334456 2.5077107686309 2.51919567722021 2.49815521640372 2.4315164539251 2.43343335040868 2.46283766412673 2.49046704878768 2.49996091899617 2.50109078577579 2.50123675502094 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000003 2.50125000000028 2.50125000000121 2.50125000014252 2.50125000862789 2.5012501270319 2.50125153878484 2.50126596357785 2.50138839362788 2.50190829852001 2.50253143276145 2.49730997484685 2.4711846540136 2.46967768037185 2.49209666704192 2.49962733410668 2.50106042062087 2.50123011664921 2.50124844112029 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000008 2.50125000000056 2.5012500000021 2.50125000033455 2.50125001327112 2.50125017370977 2.50125204497083 2.50126434758196 2.50130648551282 2.50125020444725 2.500088796806 2.49630296546738 2.49692978129782 2.50007794208536 2.50105992832656 2.50122807561851 2.50124781172452 2.50124983359067 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000001 2.5012500000002 2.50125000000134 2.50125000000873 2.50125000039353 2.50125001382068 2.50125016314226 2.50125112772164 2.50125369789792 2.50123439139653 2.50101213052607 2.50054791568097 2.50074191177358 2.50111383584011 2.50123123722341 2.50124791303752 2.50124980202498 2.50124998550472 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000001 2.50125000000025 2.50125000000158 2.50125000000948 2.501250000353 2.50125001051172 2.50125007052678 2.50125014285525 2.50124683482877 2.50122560550594 2.50116783553482 2.50119817476445 2.50123702976198 2.50124842058693 2.50124982931386 2.50124998437139 2.50124999949008 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000002 2.50125000000028 2.50125000000203 2.50125000001766 2.50125000041672 2.50125000364401 2.50124999768548 2.50124958683181 2.50124766437121 2.50124205367134 2.50124550297963 2.50124894971964 2.50124988281059 2.50124998775947 2.50124999952084 2.50124999998755 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000001 2.50125000000022 2.50125000000138 2.50125000000333 2.50125000000754 2.50124999945287 2.50124996796066 2.50124980957423 2.50124935003251 2.5012496611593 2.50124992582376 2.5012499926029 2.5012499997137 2.50124999999317 2.50124999999797 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000001 2.50125000000016 2.50125000000063 2.50124999999966 2.50124999999777 2.5012499982055 2.50124998626053 2.50124995365815 2.501249977251 2.50124999549192 2.50124999989905 2.50124999999832 2.50124999999895 2.50124999999973 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125000000001 2.50125000000006 2.50124999999983 2.50124999999906 2.5012499999884 2.50124999943726 2.50124999719567 2.50124999889295 2.50124999996683 2.50124999999928 2.50124999999959 2.50124999999984 2.50124999999998 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50124999999982 2.50124999999839 2.50124999999861 2.50124999999979 2.5012499999996 2.50124999999898 2.50124999999985 2.50124999999994 2.50124999999999 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50124999999998 2.50124999999969 2.50124999999962 2.50124999999987 2.50124999999987 2.50124999999985 2.50124999999998 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50124999999997 2.50124999999994 2.50124999999997 2.50124999999997 2.50124999999999 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 2.50125 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 +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_validator.py b/toolchain/mfc/case_validator.py index 86b1ca69cd..e2f5718a0c 100644 --- a/toolchain/mfc/case_validator.py +++ b/toolchain/mfc/case_validator.py @@ -909,6 +909,11 @@ def check_ibm(self): geometry == 2 and packing_method == 2, f"particle_cloud({i}) hemisphere-shell lattice packing is not implemented", ) + shell_axis = self.get(f"particle_cloud({i})%shell_axis", 3) + self.prohibit( + geometry == 2 and shell_axis not in [1, 2, 3], + f"particle_cloud({i})%shell_axis must be 1 (x), 2 (y), or 3 (z)", + ) if geometry == 2 and shell_outer_radius is not None and self._is_numeric(shell_outer_radius): x_centroid = self.get(f"particle_cloud({i})%x_centroid", None) y_centroid = self.get(f"particle_cloud({i})%y_centroid", None) @@ -919,32 +924,34 @@ def check_ibm(self): y_end = self.get("y_domain%end", None) z_beg = self.get("z_domain%beg", None) z_end = self.get("z_domain%end", None) - - if all(self._is_numeric(v) for v in [x_centroid, x_beg, x_end]): - self.prohibit( - x_centroid - shell_outer_radius < x_beg or x_centroid + shell_outer_radius > x_end, - f"particle_cloud({i}) hemisphere shell x-extent must lie within x_domain", - ) - if n > 0 and all(self._is_numeric(v) for v in [y_centroid, y_beg, y_end, radius]): - if p > 0: + # 2D has no z-axis; shell_axis values other than 1 (x) fall back to y, matching the + # fixed +y orientation used before shell_axis existed (see s_sample_cloud_candidate). + open_axis = shell_axis if (p > 0 or shell_axis == 1) else 2 + + axes = [ + (1, "x", x_centroid, x_beg, x_end), + (2, "y", y_centroid, y_beg, y_end), + (3, "z", z_centroid, z_beg, z_end), + ] + for axis_id, name, centroid, beg, end in axes: + if axis_id == 2 and n == 0: + continue + if axis_id == 3 and p == 0: + continue + if not all(self._is_numeric(v) for v in [centroid, beg, end, radius]): + continue + if axis_id == open_axis: + # the flat face sits at the centroid and the shell opens toward +axis; require + # one particle radius of standoff so no particle surface sits on the domain wall. self.prohibit( - y_centroid - shell_outer_radius < y_beg or y_centroid + shell_outer_radius > y_end, - f"particle_cloud({i}) hemisphere shell y-extent must lie within y_domain", + centroid - radius < beg or centroid + shell_outer_radius > end, + f"particle_cloud({i}) hemisphere shell must clear {name}_domain by one particle radius", ) else: - # 2D half-annulus opens toward +y from the flat face at y_centroid; require one - # particle radius of standoff so no particle surface sits on the domain wall. self.prohibit( - y_centroid - radius < y_beg or y_centroid + shell_outer_radius > y_end, - f"particle_cloud({i}) half-annulus must clear y_domain by one particle radius", + centroid - shell_outer_radius < beg or centroid + shell_outer_radius > end, + f"particle_cloud({i}) hemisphere shell {name}-extent must lie within {name}_domain", ) - if p > 0 and all(self._is_numeric(v) for v in [z_centroid, z_beg, z_end, radius]): - # 3D hemisphere shell opens toward +z from the flat face at z_centroid; require one - # particle radius of standoff so no particle surface sits on the domain wall. - self.prohibit( - z_centroid - radius < z_beg or z_centroid + shell_outer_radius > z_end, - f"particle_cloud({i}) hemisphere shell must clear z_domain by one particle radius", - ) num_ib_airfoils_max = get_fortran_constants().get("num_ib_airfoils_max", 5) num_stl_models_max = get_fortran_constants().get("num_stl_models_max", 10) diff --git a/toolchain/mfc/params/definitions.py b/toolchain/mfc/params/definitions.py index 1460f656d4..83c6bfa64d 100644 --- a/toolchain/mfc/params/definitions.py +++ b/toolchain/mfc/params/definitions.py @@ -1057,6 +1057,7 @@ def _load(): _pb_attrs["moving_ibm"] = (INT, _pb_tags) _pb_attrs["seed"] = (INT, _pb_tags) _pb_attrs["cloud_geometry"] = (INT, _pb_tags) + _pb_attrs["shell_axis"] = (INT, _pb_tags) _pb_attrs["packing_method"] = (INT, _pb_tags) _pb_attrs["periodic"] = (INT, _pb_tags) REGISTRY.register_family(