From b3ad7dddf05a68b2132b5c0fc29a4c6b9c94cc5d Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Sun, 3 Dec 2023 21:16:25 +0100 Subject: [PATCH 01/22] Instead of Asserting, actually include topography --- source/geometry_model/box.cc | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/source/geometry_model/box.cc b/source/geometry_model/box.cc index 68d51d76db2..c14997a5454 100644 --- a/source/geometry_model/box.cc +++ b/source/geometry_model/box.cc @@ -320,12 +320,15 @@ namespace aspect this->simulator_is_past_initialization() == false, ExcMessage("After displacement of the free surface, this function can no longer be used to determine whether a point lies in the domain or not.")); - AssertThrow(Plugins::plugin_type_matches>(this->get_initial_topography_model()), - ExcMessage("After adding topography, this function can no longer be used " - "to determine whether a point lies in the domain or not.")); + // The maximal extents of the unperturbed box domain + Point max_point = extents+box_origin; + + // Get the topography at the current point + if(Plugins::plugin_type_matches>(this->get_initial_topography_model())) + max_point = add_topography(point); for (unsigned int d = 0; d < dim; ++d) - if (point[d] > extents[d]+box_origin[d]+std::numeric_limits::epsilon()*extents[d] || + if (point[d] > max_point[d]+std::numeric_limits::epsilon()*extents[d] || point[d] < box_origin[d]-std::numeric_limits::epsilon()*extents[d]) return false; From b742de63f9457baafed18cbe9830c61d417d281a Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Sun, 3 Dec 2023 21:20:29 +0100 Subject: [PATCH 02/22] Remove Assert as topography is removed from point --- source/geometry_model/chunk.cc | 3 --- 1 file changed, 3 deletions(-) diff --git a/source/geometry_model/chunk.cc b/source/geometry_model/chunk.cc index 139e8065ba1..93ba2e35d38 100644 --- a/source/geometry_model/chunk.cc +++ b/source/geometry_model/chunk.cc @@ -701,9 +701,6 @@ namespace aspect this->simulator_is_past_initialization() == false, ExcMessage("After displacement of the free surface, this function can no longer be used to determine whether a point lies in the domain or not.")); - AssertThrow(Plugins::plugin_type_matches>(this->get_initial_topography_model()), - ExcMessage("After adding topography, this function can no longer be used to determine whether a point lies in the domain or not.")); - const Point spherical_point = manifold->pull_back(point); for (unsigned int d = 0; d < dim; ++d) From cfae329a525f214d07d09081d495216636ece870 Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Sun, 3 Dec 2023 21:23:20 +0100 Subject: [PATCH 03/22] During the first timestep, the free surface does not deform the mesh --- source/geometry_model/box.cc | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/source/geometry_model/box.cc b/source/geometry_model/box.cc index c14997a5454..0b1de7889c3 100644 --- a/source/geometry_model/box.cc +++ b/source/geometry_model/box.cc @@ -317,7 +317,8 @@ namespace aspect Box::point_is_in_domain(const Point &point) const { AssertThrow(!this->get_parameters().mesh_deformation_enabled || - this->simulator_is_past_initialization() == false, + this->simulator_is_past_initialization() == false || + this->get_timestep_number() == 0, ExcMessage("After displacement of the free surface, this function can no longer be used to determine whether a point lies in the domain or not.")); // The maximal extents of the unperturbed box domain From 882e66a9fa728341d7eda3d1ec5a93ecca078b11 Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Mon, 4 Dec 2023 10:52:43 +0100 Subject: [PATCH 04/22] Also check for initial mesh deformation, indent --- source/geometry_model/box.cc | 43 +++++++++++++++++++++++++++++++----- 1 file changed, 38 insertions(+), 5 deletions(-) diff --git a/source/geometry_model/box.cc b/source/geometry_model/box.cc index 0b1de7889c3..4d1131b507e 100644 --- a/source/geometry_model/box.cc +++ b/source/geometry_model/box.cc @@ -21,6 +21,7 @@ #include #include +#include #include #include @@ -316,16 +317,48 @@ namespace aspect bool Box::point_is_in_domain(const Point &point) const { - AssertThrow(!this->get_parameters().mesh_deformation_enabled || - this->simulator_is_past_initialization() == false || - this->get_timestep_number() == 0, - ExcMessage("After displacement of the free surface, this function can no longer be used to determine whether a point lies in the domain or not.")); + // Assert that no initial mesh deformation is applied, + // and, if time-dependent mesh deformation is enabled, + // that we are still in timestep 0, when no deformation + // has occurred yet. + const std::set mesh_deformation_boundary_ids + = this->get_mesh_deformation_handler().get_active_mesh_deformation_boundary_indicators(); + + // If mesh deformation is not allowed on any boundary, we don't have to check. + if (!mesh_deformation_boundary_ids.empty()) + { + // Get the plugins assigned to each mesh deformation boundary. + std::map> mesh_deformation_boundary_indicators_map + = this->get_mesh_deformation_handler().get_active_mesh_deformation_names(); + + // Loop over each mesh deformation boundary. + for (const types::boundary_id id : mesh_deformation_boundary_ids) + { + // Currently there is only one plugin that applies initial mesh deformation, the ascii data plugin. + // TODO also capture future initial mesh deformation plugins. + const std::vector &names = mesh_deformation_boundary_indicators_map[id]; + for (const auto &name : names) + { + if (name == "ascii data") + { + AssertThrow(false, + ExcMessage("After applying initial mesh deformation, this function can no longer be used to determine whether a point lies in the domain.")); + } + else + { + AssertThrow(this->simulator_is_past_initialization() == false || + this->get_timestep_number() == 0, + ExcMessage("After displacement of the mesh, this function can no longer be used to determine whether a point lies in the domain or not.")); + } + } + } + } // The maximal extents of the unperturbed box domain Point max_point = extents+box_origin; // Get the topography at the current point - if(Plugins::plugin_type_matches>(this->get_initial_topography_model())) + if (Plugins::plugin_type_matches>(this->get_initial_topography_model())) max_point = add_topography(point); for (unsigned int d = 0; d < dim; ++d) From 9f8508010e21d6506c98ec18e885b387dd20a772 Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Mon, 4 Dec 2023 10:53:56 +0100 Subject: [PATCH 05/22] Update comments --- source/geometry_model/box.cc | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/source/geometry_model/box.cc b/source/geometry_model/box.cc index 4d1131b507e..46ad7d1096c 100644 --- a/source/geometry_model/box.cc +++ b/source/geometry_model/box.cc @@ -354,13 +354,14 @@ namespace aspect } } - // The maximal extents of the unperturbed box domain + // The maximal extents of the unperturbed box domain. Point max_point = extents+box_origin; - // Get the topography at the current point + // Get the topography at the current point. if (Plugins::plugin_type_matches>(this->get_initial_topography_model())) max_point = add_topography(point); + // Check whether point lies within the min/max coordinates of the domain including initial topography. for (unsigned int d = 0; d < dim; ++d) if (point[d] > max_point[d]+std::numeric_limits::epsilon()*extents[d] || point[d] < box_origin[d]-std::numeric_limits::epsilon()*extents[d]) From 2b3894dab88c5dc561c2f5193d301b23c8a66248 Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Mon, 4 Dec 2023 11:31:06 +0100 Subject: [PATCH 06/22] Fix bugs --- source/geometry_model/box.cc | 24 +++++++++++++++++++----- 1 file changed, 19 insertions(+), 5 deletions(-) diff --git a/source/geometry_model/box.cc b/source/geometry_model/box.cc index 46ad7d1096c..6d01ae48c6c 100644 --- a/source/geometry_model/box.cc +++ b/source/geometry_model/box.cc @@ -321,12 +321,17 @@ namespace aspect // and, if time-dependent mesh deformation is enabled, // that we are still in timestep 0, when no deformation // has occurred yet. - const std::set mesh_deformation_boundary_ids - = this->get_mesh_deformation_handler().get_active_mesh_deformation_boundary_indicators(); // If mesh deformation is not allowed on any boundary, we don't have to check. - if (!mesh_deformation_boundary_ids.empty()) + if (this->get_parameters().mesh_deformation_enabled) { + AssertThrow(this->simulator_is_past_initialization(), + ExcMessage("Because mesh deformation is enabled, but the simulator is not past initialization, this function cannot be used to determine whether a point lies in the domain.")); + + // Get the boundaries with assigned mesh deformation. + const std::set mesh_deformation_boundary_ids + = this->get_mesh_deformation_handler().get_active_mesh_deformation_boundary_indicators(); + // Get the plugins assigned to each mesh deformation boundary. std::map> mesh_deformation_boundary_indicators_map = this->get_mesh_deformation_handler().get_active_mesh_deformation_names(); @@ -358,8 +363,17 @@ namespace aspect Point max_point = extents+box_origin; // Get the topography at the current point. - if (Plugins::plugin_type_matches>(this->get_initial_topography_model())) - max_point = add_topography(point); + if (!Plugins::plugin_type_matches>(this->get_initial_topography_model())) + { + // Get the surface x (,y) point + Point surface_point; + for (unsigned int d=0; dvalue(surface_point); + max_point[dim-1] += topo; + } // Check whether point lies within the min/max coordinates of the domain including initial topography. for (unsigned int d = 0; d < dim; ++d) From 062ef6cbfcba578c131875fb0c264c372ed21705 Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Thu, 7 Dec 2023 22:00:05 +0100 Subject: [PATCH 07/22] Swith Assert mesh def for loop over cells --- source/geometry_model/box.cc | 63 +++++++++++++++--------------------- 1 file changed, 26 insertions(+), 37 deletions(-) diff --git a/source/geometry_model/box.cc b/source/geometry_model/box.cc index 6d01ae48c6c..469baca7087 100644 --- a/source/geometry_model/box.cc +++ b/source/geometry_model/box.cc @@ -317,52 +317,41 @@ namespace aspect bool Box::point_is_in_domain(const Point &point) const { - // Assert that no initial mesh deformation is applied, - // and, if time-dependent mesh deformation is enabled, - // that we are still in timestep 0, when no deformation - // has occurred yet. - - // If mesh deformation is not allowed on any boundary, we don't have to check. + // If mesh deformation occurs over time, and we're past the + // first timestep (during which no mesh deformation is applied), + // we can no longer use the original extents of the model domain + // to check whether the point lies in the domain. + // Initial mesh deformation is applied during the first + // timestep, but this function might be called before the + // simulator is initialized and thus before we can ask whether initial + // or time-dependent mesh deformation is applied. + // Therefore, when any mesh deformation is enabled, + // we loop over all the current grid cells to see if the point + // lies within the domain. + // Do note that this function can be called before mesh deformation + // is applied in the first timestep, + // and therefore there is no guarantee that the point will lie in + // the domain after initial mesh deformation. if (this->get_parameters().mesh_deformation_enabled) { - AssertThrow(this->simulator_is_past_initialization(), - ExcMessage("Because mesh deformation is enabled, but the simulator is not past initialization, this function cannot be used to determine whether a point lies in the domain.")); - - // Get the boundaries with assigned mesh deformation. - const std::set mesh_deformation_boundary_ids - = this->get_mesh_deformation_handler().get_active_mesh_deformation_boundary_indicators(); + std::pair::active_cell_iterator, + Point> it = + GridTools::find_active_cell_around_point<> (this->get_mapping(), this->get_triangulation(), point); - // Get the plugins assigned to each mesh deformation boundary. - std::map> mesh_deformation_boundary_indicators_map - = this->get_mesh_deformation_handler().get_active_mesh_deformation_names(); + // If no cell contains the point, throw. + AssertThrow(it.first.state() == IteratorState::valid, + ExcMessage("The given point does not lie within the domain.")); - // Loop over each mesh deformation boundary. - for (const types::boundary_id id : mesh_deformation_boundary_ids) - { - // Currently there is only one plugin that applies initial mesh deformation, the ascii data plugin. - // TODO also capture future initial mesh deformation plugins. - const std::vector &names = mesh_deformation_boundary_indicators_map[id]; - for (const auto &name : names) - { - if (name == "ascii data") - { - AssertThrow(false, - ExcMessage("After applying initial mesh deformation, this function can no longer be used to determine whether a point lies in the domain.")); - } - else - { - AssertThrow(this->simulator_is_past_initialization() == false || - this->get_timestep_number() == 0, - ExcMessage("After displacement of the mesh, this function can no longer be used to determine whether a point lies in the domain or not.")); - } - } - } + // Return true, since we found a cell that contains our point. + return true; } // The maximal extents of the unperturbed box domain. Point max_point = extents+box_origin; - // Get the topography at the current point. + // If mesh deformation is not enabled, but initial topography + // was applied to the mesh, include this topography in the + // extent of the domain. if (!Plugins::plugin_type_matches>(this->get_initial_topography_model())) { // Get the surface x (,y) point From c3a3126d645f9eb5425384d1d26949e8d7998252 Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Fri, 8 Dec 2023 16:49:50 +0100 Subject: [PATCH 08/22] Only find active cell after initialization --- source/geometry_model/box.cc | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/source/geometry_model/box.cc b/source/geometry_model/box.cc index 469baca7087..40045472929 100644 --- a/source/geometry_model/box.cc +++ b/source/geometry_model/box.cc @@ -329,21 +329,25 @@ namespace aspect // we loop over all the current grid cells to see if the point // lies within the domain. // Do note that this function can be called before mesh deformation - // is applied in the first timestep, - // and therefore there is no guarantee that the point will lie in - // the domain after initial mesh deformation. - if (this->get_parameters().mesh_deformation_enabled) + // is applied in the first timestep, and therefore there is no guarantee + // that the point will lie in the domain after initial mesh deformation. + // For example, boundary traction plugins are initialized before + // mesh deformation plugins and the initial_lithostatic_pressure + // plugin calls this function while the simulator is being initialized. + // However, if the simulator has not been initialized, + // no mesh deformation has been applied and we can use the geometry model's + // extents to check whether a point lies in the domain. + if (this->get_parameters().mesh_deformation_enabled && + this->simulator_is_past_initialization()) { std::pair::active_cell_iterator, Point> it = GridTools::find_active_cell_around_point<> (this->get_mapping(), this->get_triangulation(), point); - // If no cell contains the point, throw. - AssertThrow(it.first.state() == IteratorState::valid, - ExcMessage("The given point does not lie within the domain.")); - - // Return true, since we found a cell that contains our point. - return true; + if (it.first.state() == IteratorState::valid) + return true; + else + return false; } // The maximal extents of the unperturbed box domain. From ab980777b6047e3895f2b7f88b7e4138a8bea0c1 Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Fri, 8 Dec 2023 16:50:09 +0100 Subject: [PATCH 09/22] Also find active cell for chunk --- source/geometry_model/chunk.cc | 23 +++++++++++++++++++---- 1 file changed, 19 insertions(+), 4 deletions(-) diff --git a/source/geometry_model/chunk.cc b/source/geometry_model/chunk.cc index 93ba2e35d38..4fb70dec3eb 100644 --- a/source/geometry_model/chunk.cc +++ b/source/geometry_model/chunk.cc @@ -697,10 +697,25 @@ namespace aspect bool Chunk::point_is_in_domain(const Point &point) const { - AssertThrow(!this->get_parameters().mesh_deformation_enabled || - this->simulator_is_past_initialization() == false, - ExcMessage("After displacement of the free surface, this function can no longer be used to determine whether a point lies in the domain or not.")); - + // If mesh deformation has possibly been applied, + // e.g. initial mesh deformation after the simulator + // has been initialized or time-dependent mesh deformation + // after timestep 0, we cannot use information from the + // geometry model to determine whether the point lies in + // the model domain. Instead, we loop over all cells + // to check whether the point lies in one of them. + if (this->get_parameters().mesh_deformation_enabled && + this->simulator_is_past_initialization()) + { + std::pair::active_cell_iterator, + Point> it = + GridTools::find_active_cell_around_point<> (this->get_mapping(), this->get_triangulation(), point); + + if (it.first.state() == IteratorState::valid) + return true; + else + return false; + } const Point spherical_point = manifold->pull_back(point); for (unsigned int d = 0; d < dim; ++d) From 4c8ebf386734696a0c4ab95f5fd4bacb98558bc4 Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Fri, 8 Dec 2023 17:19:34 +0100 Subject: [PATCH 10/22] Add test box initial topo mesh def --- ...tasy_initial_topo_lithostatic_pressure.prm | 40 ++++++++++++++++ .../screen-output | 46 +++++++++++++++++++ .../statistics | 20 ++++++++ 3 files changed, 106 insertions(+) create mode 100644 tests/airy_isostasy_initial_topo_lithostatic_pressure.prm create mode 100644 tests/airy_isostasy_initial_topo_lithostatic_pressure/screen-output create mode 100644 tests/airy_isostasy_initial_topo_lithostatic_pressure/statistics diff --git a/tests/airy_isostasy_initial_topo_lithostatic_pressure.prm b/tests/airy_isostasy_initial_topo_lithostatic_pressure.prm new file mode 100644 index 00000000000..c690e7e9e24 --- /dev/null +++ b/tests/airy_isostasy_initial_topo_lithostatic_pressure.prm @@ -0,0 +1,40 @@ +# We setup three columns, of which the middle column has a higher density +# and viscosity. By choosing a free surface in combination with a +# prescribed lithostatic pressure, +# we see the middle column sink until isostasy is reached (if end time is increased) +# Compared to airy_isostasy.prm, we start off with a nonzero initial topography. +# This test in particular tests the point_is_in_domain function of the geometry +# model. It should take into account the initial topography and find that the +# Representative point lies within the domain. At x = 0.2, the initial topography +# is namely 0.02948. + +include $ASPECT_SOURCE_DIR/tests/airy_isostasy.prm + +set Dimension = 2 +# The point_is_in_domain function is only called in +# timestep 0, so we do not need to run the test further. +set End time = 0 + +subsection Geometry model + subsection Initial topography model + set Model name = function + + subsection Function + set Function constants = L=1, R=0.05, C=0.5 + set Function expression = R * cos(pi*(x-C)/(L)) + set Maximum topography value = 0.05 + end + end +end + +# The representative point determines where the pressure +# profile is calculated. We pick a point in the left +# column which has some initial topography prescribed. +subsection Boundary traction model + set Prescribed traction boundary indicators = 2 y: initial lithostatic pressure + + subsection Initial lithostatic pressure + set Number of integration points = 260 + set Representative point = 0.2,1.02 + end +end diff --git a/tests/airy_isostasy_initial_topo_lithostatic_pressure/screen-output b/tests/airy_isostasy_initial_topo_lithostatic_pressure/screen-output new file mode 100644 index 00000000000..e5d91664954 --- /dev/null +++ b/tests/airy_isostasy_initial_topo_lithostatic_pressure/screen-output @@ -0,0 +1,46 @@ + +Number of active cells: 1,024 (on 6 levels) +Number of degrees of freedom: 17,989 (8,450+1,089+4,225+4,225) + +Number of mesh deformation degrees of freedom: 2,178 + Solving mesh displacement system... 0 iterations. +*** Timestep 0: t=0 years, dt=0 years + Solving mesh displacement system... 0 iterations. + Skipping temperature solve because RHS is zero. + Solving C_1 system ... 0 iterations. + Rebuilding Stokes preconditioner... + Solving Stokes system... 200+8 iterations. + +Number of active cells: 1,204 (on 7 levels) +Number of degrees of freedom: 21,360 (10,034+1,292+5,017+5,017) + +Number of mesh deformation degrees of freedom: 2,584 + Solving mesh displacement system... 0 iterations. +*** Timestep 0: t=0 years, dt=0 years + Solving mesh displacement system... 0 iterations. + Skipping temperature solve because RHS is zero. + Solving C_1 system ... 0 iterations. + Rebuilding Stokes preconditioner... + Solving Stokes system... 167+0 iterations. + +Number of active cells: 1,948 (on 8 levels) +Number of degrees of freedom: 34,621 (16,270+2,081+8,135+8,135) + +Number of mesh deformation degrees of freedom: 4,162 + Solving mesh displacement system... 0 iterations. +*** Timestep 0: t=0 years, dt=0 years + Solving mesh displacement system... 0 iterations. + Skipping temperature solve because RHS is zero. + Solving C_1 system ... 0 iterations. + Rebuilding Stokes preconditioner... + Solving Stokes system... 200+6 iterations. + + Postprocessing: + RMS, max velocity: 6.63e+05 m/year, 8.88e+05 m/year + Pressure min/avg/max: -0.6507 Pa, 5.005 Pa, 10 Pa + Topography min/max: 0 m, 0 m + +Termination requested by criterion: end time + + + diff --git a/tests/airy_isostasy_initial_topo_lithostatic_pressure/statistics b/tests/airy_isostasy_initial_topo_lithostatic_pressure/statistics new file mode 100644 index 00000000000..70cc69bf902 --- /dev/null +++ b/tests/airy_isostasy_initial_topo_lithostatic_pressure/statistics @@ -0,0 +1,20 @@ +# 1: Time step number +# 2: Time (years) +# 3: Time step size (years) +# 4: Number of mesh cells +# 5: Number of Stokes degrees of freedom +# 6: Number of temperature degrees of freedom +# 7: Number of degrees of freedom for all compositions +# 8: Iterations for temperature solver +# 9: Iterations for composition solver 1 +# 10: Iterations for Stokes solver +# 11: Velocity iterations in Stokes preconditioner +# 12: Schur complement iterations in Stokes preconditioner +# 13: RMS velocity (m/year) +# 14: Max. velocity (m/year) +# 15: Minimal pressure (Pa) +# 16: Average pressure (Pa) +# 17: Maximal pressure (Pa) +# 18: Minimum topography (m) +# 19: Maximum topography (m) +0 0.000000000000e+00 0.000000000000e+00 1948 18351 8135 8135 0 0 206 444 840 6.62545995e+05 8.87572504e+05 -6.50724076e-01 5.00501899e+00 1.00000000e+01 0.00000000e+00 0.00000000e+00 From e17370a330deac165e99bf7a1b02abb01f6f8b7e Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Fri, 8 Dec 2023 17:33:09 +0100 Subject: [PATCH 11/22] Add test chunk initial topo lithostatic pressure --- ...ostatic_pressure_2d_initial_topography.prm | 45 ++ .../screen-output | 19 + .../solution/solution-00000.0000.gnuplot | 391 ++++++++++++++++++ .../statistics | 19 + 4 files changed, 474 insertions(+) create mode 100644 tests/chunk_lithostatic_pressure_2d_initial_topography.prm create mode 100644 tests/chunk_lithostatic_pressure_2d_initial_topography/screen-output create mode 100644 tests/chunk_lithostatic_pressure_2d_initial_topography/solution/solution-00000.0000.gnuplot create mode 100644 tests/chunk_lithostatic_pressure_2d_initial_topography/statistics diff --git a/tests/chunk_lithostatic_pressure_2d_initial_topography.prm b/tests/chunk_lithostatic_pressure_2d_initial_topography.prm new file mode 100644 index 00000000000..7861c6a79ae --- /dev/null +++ b/tests/chunk_lithostatic_pressure_2d_initial_topography.prm @@ -0,0 +1,45 @@ +# In this variation on chunk_lithostatic_pressure_2d.prm +# we test whether the point_is_in_domain function +# handles the initial topography correctly. We therefore +# add an initial topography and pick a representative point +# for the initial lithostatic pressure that lies above +# the unperturbed surface but below the surface plus +# initial topography. + +include $ASPECT_SOURCE_DIR/tests/chunk_lithostatic_pressure_2d.prm + +set Dimension = 3 +set End time = 0 + +subsection Geometry model + set Model name = chunk + + subsection Chunk + set Chunk inner radius = 3471000 + set Chunk outer radius = 6371000 + set Chunk minimum longitude = 10 + set Chunk maximum longitude = 50 + set Chunk minimum latitude = 0 + set Chunk maximum latitude = 40 + end + + subsection Initial topography model + set Model name = function + + subsection Function + set Function expression = 5000 + set Maximum topography value = 5000 + end + end +end + + +subsection Boundary traction model + set Prescribed traction boundary indicators = east: initial lithostatic pressure + + # Set a radius between 6371 and 6376 km (surface radius + # without and with initial topography. + subsection Initial lithostatic pressure + set Representative point = 6374000, 50, 40 + end +end diff --git a/tests/chunk_lithostatic_pressure_2d_initial_topography/screen-output b/tests/chunk_lithostatic_pressure_2d_initial_topography/screen-output new file mode 100644 index 00000000000..3fc22aca562 --- /dev/null +++ b/tests/chunk_lithostatic_pressure_2d_initial_topography/screen-output @@ -0,0 +1,19 @@ + +Number of active cells: 8 (on 2 levels) +Number of degrees of freedom: 652 (375+27+125+125) + +*** Timestep 0: t=0 years, dt=0 years + Solving temperature system... 0 iterations. + Solving C_1 system ... 11 iterations. + Rebuilding Stokes preconditioner... + Solving Stokes system... 22+0 iterations. + + Postprocessing: + RMS, max velocity: 0.309 m/year, 1.86 m/year + Pressure min/avg/max: -1.965e+07 Pa, 3.29e+10 Pa, 9.157e+10 Pa + Writing graphical output: output-chunk_lithostatic_pressure_2d/solution/solution-00000 + +Termination requested by criterion: end time + + + diff --git a/tests/chunk_lithostatic_pressure_2d_initial_topography/solution/solution-00000.0000.gnuplot b/tests/chunk_lithostatic_pressure_2d_initial_topography/solution/solution-00000.0000.gnuplot new file mode 100644 index 00000000000..baca9c2fcbc --- /dev/null +++ b/tests/chunk_lithostatic_pressure_2d_initial_topography/solution/solution-00000.0000.gnuplot @@ -0,0 +1,391 @@ +# This file was generated by the deal.II library. + + +# +# For a description of the GNUPLOT format see the GNUPLOT manual. +# +#

+3.41827e+06 602733 0 0 0 0 9.14242e+10 1600 0 9.14598e+10 4.13146e+06 3.14744e+07 4.13146e+06 9.14342e+10 5.54979e+06 3.14744e+07 5.54979e+06 9.10905e+10 1.04362e-13 4.50663e+07 +4.84624e+06 854523 0 -0.057618 -0.0101596 0 4.56474e+10 1600 -1.84884e-16 4.56173e+10 -6.04158e+06 8.6398e+06 -6.04158e+06 4.56496e+10 1.52343e+06 8.6398e+06 1.52343e+06 4.56317e+10 9.68689e-15 -1.78755e+07 + + +3.41827e+06 602733 0 0 0 0 9.14242e+10 1600 0 9.14598e+10 4.13146e+06 3.14744e+07 4.13146e+06 9.14342e+10 5.54979e+06 3.14744e+07 5.54979e+06 9.10905e+10 1.04362e-13 4.50663e+07 +3.00597e+06 1.7355e+06 0 0.0355464 -0.0615682 0 9.14379e+10 1600 0 9.1461e+10 1.6499e+07 2.71656e+07 1.6499e+07 9.1447e+10 1.7322e+07 2.71656e+07 1.7322e+07 9.1101e+10 1.03565e-13 5.8845e+07 + + +3.41827e+06 602733 0 0 0 0 9.14242e+10 1600 0 9.14598e+10 4.13146e+06 3.14744e+07 4.13146e+06 9.14342e+10 5.54979e+06 3.14744e+07 5.54979e+06 9.10905e+10 1.04362e-13 4.50663e+07 +3.21212e+06 566384 1.18715e+06 0.185945 0.0327872 -0.518762 9.13123e+10 1600 0 9.14457e+10 -3.67688e+07 -1.49744e+08 -3.67688e+07 9.12698e+10 1.57118e+08 -1.49744e+08 1.57118e+08 9.16547e+10 1.46277e-13 -6.68253e+07 + + +4.84624e+06 854523 0 -0.057618 -0.0101596 0 4.56474e+10 1600 -1.84884e-16 4.56173e+10 -6.04158e+06 8.6398e+06 -6.04158e+06 4.56496e+10 1.52343e+06 8.6398e+06 1.52343e+06 4.56317e+10 9.68689e-15 -1.78755e+07 +4.26171e+06 2.4605e+06 0 -0.000689281 -0.112987 0 4.56467e+10 1600 -9.13492e-17 4.56287e+10 -1.6118e+07 7.29339e+06 -1.6118e+07 4.56415e+10 4.68681e+06 7.29339e+06 4.68681e+06 4.56313e+10 9.75342e-15 -1.85507e+07 + + +4.84624e+06 854523 0 -0.057618 -0.0101596 0 4.56474e+10 1600 -1.84884e-16 4.56173e+10 -6.04158e+06 8.6398e+06 -6.04158e+06 4.56496e+10 1.52343e+06 8.6398e+06 1.52343e+06 4.56317e+10 9.68689e-15 -1.78755e+07 +4.55397e+06 802989 1.68308e+06 -0.148579 -0.0261985 -0.247956 4.56703e+10 1600 -5.52646e-17 4.56203e+10 -1.9544e+07 -2.35491e+06 -1.9544e+07 4.56633e+10 -223785 -2.35491e+06 -223785 4.57163e+10 2.59646e-14 5.05121e+06 + + +3.00597e+06 1.7355e+06 0 0.0355464 -0.0615682 0 9.14379e+10 1600 0 9.1461e+10 1.6499e+07 2.71656e+07 1.6499e+07 9.1447e+10 1.7322e+07 2.71656e+07 1.7322e+07 9.1101e+10 1.03565e-13 5.8845e+07 +4.26171e+06 2.4605e+06 0 -0.000689281 -0.112987 0 4.56467e+10 1600 -9.13492e-17 4.56287e+10 -1.6118e+07 7.29339e+06 -1.6118e+07 4.56415e+10 4.68681e+06 7.29339e+06 4.68681e+06 4.56313e+10 9.75342e-15 -1.85507e+07 + + +3.00597e+06 1.7355e+06 0 0.0355464 -0.0615682 0 9.14379e+10 1600 0 9.1461e+10 1.6499e+07 2.71656e+07 1.6499e+07 9.1447e+10 1.7322e+07 2.71656e+07 1.7322e+07 9.1101e+10 1.03565e-13 5.8845e+07 +2.82469e+06 1.63084e+06 1.18715e+06 0.19034 0.0376886 -0.504666 9.13093e+10 1600 0 9.13328e+10 9.30997e+07 -1.30653e+07 9.30997e+07 9.13805e+10 -2.1834e+08 -1.30653e+07 -2.1834e+08 9.16534e+10 1.46993e-13 -6.9817e+07 + + +4.26171e+06 2.4605e+06 0 -0.000689281 -0.112987 0 4.56467e+10 1600 -9.13492e-17 4.56287e+10 -1.6118e+07 7.29339e+06 -1.6118e+07 4.56415e+10 4.68681e+06 7.29339e+06 4.68681e+06 4.56313e+10 9.75342e-15 -1.85507e+07 +4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56204e+10 -2.00601e+07 684454 -2.00601e+07 4.56656e+10 -4.28112e+06 684454 -4.28112e+06 4.57146e+10 2.56922e-14 4.8163e+06 + + +3.21212e+06 566384 1.18715e+06 0.185945 0.0327872 -0.518762 9.13123e+10 1600 0 9.14457e+10 -3.67688e+07 -1.49744e+08 -3.67688e+07 9.12698e+10 1.57118e+08 -1.49744e+08 1.57118e+08 9.16547e+10 1.46277e-13 -6.68253e+07 +4.55397e+06 802989 1.68308e+06 -0.148579 -0.0261985 -0.247956 4.56703e+10 1600 -5.52646e-17 4.56203e+10 -1.9544e+07 -2.35491e+06 -1.9544e+07 4.56633e+10 -223785 -2.35491e+06 -223785 4.57163e+10 2.59646e-14 5.05121e+06 + + +3.21212e+06 566384 1.18715e+06 0.185945 0.0327872 -0.518762 9.13123e+10 1600 0 9.14457e+10 -3.67688e+07 -1.49744e+08 -3.67688e+07 9.12698e+10 1.57118e+08 -1.49744e+08 1.57118e+08 9.16547e+10 1.46277e-13 -6.68253e+07 +2.82469e+06 1.63084e+06 1.18715e+06 0.19034 0.0376886 -0.504666 9.13093e+10 1600 0 9.13328e+10 9.30997e+07 -1.30653e+07 9.30997e+07 9.13805e+10 -2.1834e+08 -1.30653e+07 -2.1834e+08 9.16534e+10 1.46993e-13 -6.9817e+07 + + +4.55397e+06 802989 1.68308e+06 -0.148579 -0.0261985 -0.247956 4.56703e+10 1600 -5.52646e-17 4.56203e+10 -1.9544e+07 -2.35491e+06 -1.9544e+07 4.56633e+10 -223785 -2.35491e+06 -223785 4.57163e+10 2.59646e-14 5.05121e+06 +4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56204e+10 -2.00601e+07 684454 -2.00601e+07 4.56656e+10 -4.28112e+06 684454 -4.28112e+06 4.57146e+10 2.56922e-14 4.8163e+06 + + +2.82469e+06 1.63084e+06 1.18715e+06 0.19034 0.0376886 -0.504666 9.13093e+10 1600 0 9.13328e+10 9.30997e+07 -1.30653e+07 9.30997e+07 9.13805e+10 -2.1834e+08 -1.30653e+07 -2.1834e+08 9.16534e+10 1.46993e-13 -6.9817e+07 +4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56204e+10 -2.00601e+07 684454 -2.00601e+07 4.56656e+10 -4.28112e+06 684454 -4.28112e+06 4.57146e+10 2.56922e-14 4.8163e+06 + + +4.84624e+06 854523 0 -0.057618 -0.0101596 0 4.56474e+10 1600 -1.84884e-16 4.56422e+10 -1.65431e+06 8.63978e+06 -1.65431e+06 4.56504e+10 1.52343e+06 8.63978e+06 1.52343e+06 4.56317e+10 6.45996e-15 -1.78755e+07 +6.27421e+06 1.10631e+06 0 0 0 0 6.38156e+06 1600 1 4.80547e+08 -1.17027e+08 9.8169e+07 -1.17027e+08 6.31248e+08 1.73098e+07 9.8169e+07 1.73098e+07 3.1943e+08 1.09486e-16 6.38156e+06 + + +4.84624e+06 854523 0 -0.057618 -0.0101596 0 4.56474e+10 1600 -1.84884e-16 4.56422e+10 -1.65431e+06 8.63978e+06 -1.65431e+06 4.56504e+10 1.52343e+06 8.63978e+06 1.52343e+06 4.56317e+10 6.45996e-15 -1.78755e+07 +4.26171e+06 2.4605e+06 0 -0.000689281 -0.112987 0 4.56467e+10 1600 -9.13492e-17 4.56486e+10 -6.88188e+06 7.29337e+06 -6.88188e+06 4.56455e+10 4.68681e+06 7.29337e+06 4.68681e+06 4.56313e+10 7.1962e-15 -1.85507e+07 + + +4.84624e+06 854523 0 -0.057618 -0.0101596 0 4.56474e+10 1600 -1.84884e-16 4.56422e+10 -1.65431e+06 8.63978e+06 -1.65431e+06 4.56504e+10 1.52343e+06 8.63978e+06 1.52343e+06 4.56317e+10 6.45996e-15 -1.78755e+07 +4.55397e+06 802989 1.68308e+06 -0.148579 -0.0261985 -0.247956 4.56703e+10 1600 -5.52646e-17 4.56596e+10 -1.26263e+07 -2.35403e+07 -1.26263e+07 4.56645e+10 -3.95929e+06 -2.35403e+07 -3.95929e+06 4.56952e+10 1.66034e-14 5.05121e+06 + + +6.27421e+06 1.10631e+06 0 0 0 0 6.38156e+06 1600 1 4.80547e+08 -1.17027e+08 9.8169e+07 -1.17027e+08 6.31248e+08 1.73098e+07 9.8169e+07 1.73098e+07 3.1943e+08 1.09486e-16 6.38156e+06 +5.51745e+06 3.1855e+06 0 -0.00516999 0.00895468 0 1.11727e+07 1600 1 -1.34887e+08 3.51473e+08 6.55746e+07 3.51473e+08 -1.22425e+09 5.55838e+07 6.55746e+07 5.55838e+07 2.56716e+08 4.24209e-16 1.11727e+07 + + +6.27421e+06 1.10631e+06 0 0 0 0 6.38156e+06 1600 1 4.80547e+08 -1.17027e+08 9.8169e+07 -1.17027e+08 6.31248e+08 1.73098e+07 9.8169e+07 1.73098e+07 3.1943e+08 1.09486e-16 6.38156e+06 +5.89583e+06 1.03959e+06 2.17901e+06 0.0170632 0.00300871 -0.0476041 -3.4701e+06 1600 1 5.00929e+09 1.10033e+09 1.30345e+09 1.10033e+09 -1.16799e+09 1.40563e+08 1.30345e+09 1.40563e+08 3.27351e+09 1.80826e-15 -3.4701e+06 + + +4.26171e+06 2.4605e+06 0 -0.000689281 -0.112987 0 4.56467e+10 1600 -9.13492e-17 4.56486e+10 -6.88188e+06 7.29337e+06 -6.88188e+06 4.56455e+10 4.68681e+06 7.29337e+06 4.68681e+06 4.56313e+10 7.1962e-15 -1.85507e+07 +5.51745e+06 3.1855e+06 0 -0.00516999 0.00895468 0 1.11727e+07 1600 1 -1.34887e+08 3.51473e+08 6.55746e+07 3.51473e+08 -1.22425e+09 5.55838e+07 6.55746e+07 5.55838e+07 2.56716e+08 4.24209e-16 1.11727e+07 + + +4.26171e+06 2.4605e+06 0 -0.000689281 -0.112987 0 4.56467e+10 1600 -9.13492e-17 4.56486e+10 -6.88188e+06 7.29337e+06 -6.88188e+06 4.56455e+10 4.68681e+06 7.29337e+06 4.68681e+06 4.56313e+10 7.1962e-15 -1.85507e+07 +4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56509e+10 -3.6295e+06 -1.76123e+07 -3.6295e+06 4.56744e+10 -1.53349e+07 -1.76123e+07 -1.53349e+07 4.56938e+10 1.59777e-14 4.8163e+06 + + +5.51745e+06 3.1855e+06 0 -0.00516999 0.00895468 0 1.11727e+07 1600 1 -1.34887e+08 3.51473e+08 6.55746e+07 3.51473e+08 -1.22425e+09 5.55838e+07 6.55746e+07 5.55838e+07 2.56716e+08 4.24209e-16 1.11727e+07 +5.18471e+06 2.99339e+06 2.17901e+06 -0.0171684 0.0586576 -0.0397298 -1.83017e+06 1600 1 3.17549e+09 3.75887e+09 2.18939e+09 3.75887e+09 -3.25586e+08 -6.23188e+08 2.18939e+09 -6.23188e+08 2.94865e+09 2.40568e-15 -1.83017e+06 + + +4.55397e+06 802989 1.68308e+06 -0.148579 -0.0261985 -0.247956 4.56703e+10 1600 -5.52646e-17 4.56596e+10 -1.26263e+07 -2.35403e+07 -1.26263e+07 4.56645e+10 -3.95929e+06 -2.35403e+07 -3.95929e+06 4.56952e+10 1.66034e-14 5.05121e+06 +5.89583e+06 1.03959e+06 2.17901e+06 0.0170632 0.00300871 -0.0476041 -3.4701e+06 1600 1 5.00929e+09 1.10033e+09 1.30345e+09 1.10033e+09 -1.16799e+09 1.40563e+08 1.30345e+09 1.40563e+08 3.27351e+09 1.80826e-15 -3.4701e+06 + + +4.55397e+06 802989 1.68308e+06 -0.148579 -0.0261985 -0.247956 4.56703e+10 1600 -5.52646e-17 4.56596e+10 -1.26263e+07 -2.35403e+07 -1.26263e+07 4.56645e+10 -3.95929e+06 -2.35403e+07 -3.95929e+06 4.56952e+10 1.66034e-14 5.05121e+06 +4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56509e+10 -3.6295e+06 -1.76123e+07 -3.6295e+06 4.56744e+10 -1.53349e+07 -1.76123e+07 -1.53349e+07 4.56938e+10 1.59777e-14 4.8163e+06 + + +5.89583e+06 1.03959e+06 2.17901e+06 0.0170632 0.00300871 -0.0476041 -3.4701e+06 1600 1 5.00929e+09 1.10033e+09 1.30345e+09 1.10033e+09 -1.16799e+09 1.40563e+08 1.30345e+09 1.40563e+08 3.27351e+09 1.80826e-15 -3.4701e+06 +5.18471e+06 2.99339e+06 2.17901e+06 -0.0171684 0.0586576 -0.0397298 -1.83017e+06 1600 1 3.17549e+09 3.75887e+09 2.18939e+09 3.75887e+09 -3.25586e+08 -6.23188e+08 2.18939e+09 -6.23188e+08 2.94865e+09 2.40568e-15 -1.83017e+06 + + +4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56509e+10 -3.6295e+06 -1.76123e+07 -3.6295e+06 4.56744e+10 -1.53349e+07 -1.76123e+07 -1.53349e+07 4.56938e+10 1.59777e-14 4.8163e+06 +5.18471e+06 2.99339e+06 2.17901e+06 -0.0171684 0.0586576 -0.0397298 -1.83017e+06 1600 1 3.17549e+09 3.75887e+09 2.18939e+09 3.75887e+09 -3.25586e+08 -6.23188e+08 2.18939e+09 -6.23188e+08 2.94865e+09 2.40568e-15 -1.83017e+06 + + +3.00597e+06 1.7355e+06 0 0.0355464 -0.0615682 0 9.14379e+10 1600 0 9.1468e+10 7.93689e+06 2.71656e+07 7.93689e+06 9.14555e+10 1.7322e+07 2.71656e+07 1.7322e+07 9.1101e+10 1.05519e-13 5.8845e+07 +4.26171e+06 2.4605e+06 0 -0.000689281 -0.112987 0 4.56467e+10 1600 -9.13492e-17 4.56324e+10 -1.8969e+07 7.29339e+06 -1.8969e+07 4.564e+10 4.68681e+06 7.29339e+06 4.68681e+06 4.56313e+10 1.06925e-14 -1.85507e+07 + + +3.00597e+06 1.7355e+06 0 0.0355464 -0.0615682 0 9.14379e+10 1600 0 9.1468e+10 7.93689e+06 2.71656e+07 7.93689e+06 9.14555e+10 1.7322e+07 2.71656e+07 1.7322e+07 9.1101e+10 1.05519e-13 5.8845e+07 +2.23112e+06 2.65894e+06 0 -0.425883 0.357358 0 9.13978e+10 1600 0 9.13185e+10 5.12101e+07 761904 5.12101e+07 9.14355e+10 4.14712e+07 761904 4.14712e+07 9.10642e+10 1.00487e-13 1.87398e+07 + + +3.00597e+06 1.7355e+06 0 0.0355464 -0.0615682 0 9.14379e+10 1600 0 9.1468e+10 7.93689e+06 2.71656e+07 7.93689e+06 9.14555e+10 1.7322e+07 2.71656e+07 1.7322e+07 9.1101e+10 1.05519e-13 5.8845e+07 +2.82469e+06 1.63084e+06 1.18715e+06 0.19034 0.0376886 -0.504666 9.13093e+10 1600 0 9.14467e+10 2.54924e+07 -1.92175e+08 2.54924e+07 9.12729e+10 9.18863e+07 -1.92175e+08 9.18863e+07 9.16534e+10 1.43442e-13 -6.9817e+07 + + +4.26171e+06 2.4605e+06 0 -0.000689281 -0.112987 0 4.56467e+10 1600 -9.13492e-17 4.56324e+10 -1.8969e+07 7.29339e+06 -1.8969e+07 4.564e+10 4.68681e+06 7.29339e+06 4.68681e+06 4.56313e+10 1.06925e-14 -1.85507e+07 +3.16316e+06 3.7697e+06 0 0.0866503 -0.288584 0 4.56761e+10 1600 0 4.56746e+10 -1.27688e+07 5.09732e+06 -1.27688e+07 4.56554e+10 6.93062e+06 5.09732e+06 6.93062e+06 4.56659e+10 9.07151e-15 1.08148e+07 + + +4.26171e+06 2.4605e+06 0 -0.000689281 -0.112987 0 4.56467e+10 1600 -9.13492e-17 4.56324e+10 -1.8969e+07 7.29339e+06 -1.8969e+07 4.564e+10 4.68681e+06 7.29339e+06 4.68681e+06 4.56313e+10 1.06925e-14 -1.85507e+07 +4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56466e+10 -3.70481e+07 989103 -3.70481e+07 4.56457e+10 -4.80879e+06 989103 -4.80879e+06 4.57146e+10 2.71892e-14 4.8163e+06 + + +2.23112e+06 2.65894e+06 0 -0.425883 0.357358 0 9.13978e+10 1600 0 9.13185e+10 5.12101e+07 761904 5.12101e+07 9.14355e+10 4.14712e+07 761904 4.14712e+07 9.10642e+10 1.00487e-13 1.87398e+07 +3.16316e+06 3.7697e+06 0 0.0866503 -0.288584 0 4.56761e+10 1600 0 4.56746e+10 -1.27688e+07 5.09732e+06 -1.27688e+07 4.56554e+10 6.93062e+06 5.09732e+06 6.93062e+06 4.56659e+10 9.07151e-15 1.08148e+07 + + +2.23112e+06 2.65894e+06 0 -0.425883 0.357358 0 9.13978e+10 1600 0 9.13185e+10 5.12101e+07 761904 5.12101e+07 9.14355e+10 4.14712e+07 761904 4.14712e+07 9.10642e+10 1.00487e-13 1.87398e+07 +2.09656e+06 2.49859e+06 1.18715e+06 0.151969 0.147411 -0.578639 9.13366e+10 1600 0 9.12783e+10 6.4997e+07 4.33178e+07 6.4997e+07 9.14735e+10 -1.96226e+08 4.33178e+07 -1.96226e+08 9.16845e+10 1.46516e-13 -4.25455e+07 + + +3.16316e+06 3.7697e+06 0 0.0866503 -0.288584 0 4.56761e+10 1600 0 4.56746e+10 -1.27688e+07 5.09732e+06 -1.27688e+07 4.56554e+10 6.93062e+06 5.09732e+06 6.93062e+06 4.56659e+10 9.07151e-15 1.08148e+07 +2.9724e+06 3.54236e+06 1.68308e+06 -0.120017 -0.203298 -0.283249 4.56924e+10 1600 0 4.56678e+10 -1.6427e+07 8.96147e+06 -1.6427e+07 4.56519e+10 -1.03223e+07 8.96147e+06 -1.03223e+07 4.57353e+10 2.45702e-14 2.71349e+07 + + +2.82469e+06 1.63084e+06 1.18715e+06 0.19034 0.0376886 -0.504666 9.13093e+10 1600 0 9.14467e+10 2.54924e+07 -1.92175e+08 2.54924e+07 9.12729e+10 9.18863e+07 -1.92175e+08 9.18863e+07 9.16534e+10 1.43442e-13 -6.9817e+07 +4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56466e+10 -3.70481e+07 989103 -3.70481e+07 4.56457e+10 -4.80879e+06 989103 -4.80879e+06 4.57146e+10 2.71892e-14 4.8163e+06 + + +2.82469e+06 1.63084e+06 1.18715e+06 0.19034 0.0376886 -0.504666 9.13093e+10 1600 0 9.14467e+10 2.54924e+07 -1.92175e+08 2.54924e+07 9.12729e+10 9.18863e+07 -1.92175e+08 9.18863e+07 9.16534e+10 1.43442e-13 -6.9817e+07 +2.09656e+06 2.49859e+06 1.18715e+06 0.151969 0.147411 -0.578639 9.13366e+10 1600 0 9.12783e+10 6.4997e+07 4.33178e+07 6.4997e+07 9.14735e+10 -1.96226e+08 4.33178e+07 -1.96226e+08 9.16845e+10 1.46516e-13 -4.25455e+07 + + +4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56466e+10 -3.70481e+07 989103 -3.70481e+07 4.56457e+10 -4.80879e+06 989103 -4.80879e+06 4.57146e+10 2.71892e-14 4.8163e+06 +2.9724e+06 3.54236e+06 1.68308e+06 -0.120017 -0.203298 -0.283249 4.56924e+10 1600 0 4.56678e+10 -1.6427e+07 8.96147e+06 -1.6427e+07 4.56519e+10 -1.03223e+07 8.96147e+06 -1.03223e+07 4.57353e+10 2.45702e-14 2.71349e+07 + + +2.09656e+06 2.49859e+06 1.18715e+06 0.151969 0.147411 -0.578639 9.13366e+10 1600 0 9.12783e+10 6.4997e+07 4.33178e+07 6.4997e+07 9.14735e+10 -1.96226e+08 4.33178e+07 -1.96226e+08 9.16845e+10 1.46516e-13 -4.25455e+07 +2.9724e+06 3.54236e+06 1.68308e+06 -0.120017 -0.203298 -0.283249 4.56924e+10 1600 0 4.56678e+10 -1.6427e+07 8.96147e+06 -1.6427e+07 4.56519e+10 -1.03223e+07 8.96147e+06 -1.03223e+07 4.57353e+10 2.45702e-14 2.71349e+07 + + +4.26171e+06 2.4605e+06 0 -0.000689281 -0.112987 0 4.56467e+10 1600 -9.13492e-17 4.56524e+10 -9.73295e+06 7.29337e+06 -9.73295e+06 4.5644e+10 4.68681e+06 7.29337e+06 4.68681e+06 4.56313e+10 8.39786e-15 -1.85507e+07 +5.51745e+06 3.1855e+06 0 -0.00516999 0.00895468 0 1.11727e+07 1600 1 -3.23725e+07 1.88972e+08 6.55746e+07 1.88972e+08 -9.6888e+08 5.55838e+07 6.55746e+07 5.55838e+07 2.56716e+08 3.36731e-16 1.11727e+07 + + +4.26171e+06 2.4605e+06 0 -0.000689281 -0.112987 0 4.56467e+10 1600 -9.13492e-17 4.56524e+10 -9.73295e+06 7.29337e+06 -9.73295e+06 4.5644e+10 4.68681e+06 7.29337e+06 4.68681e+06 4.56313e+10 8.39786e-15 -1.85507e+07 +3.16316e+06 3.7697e+06 0 0.0866503 -0.288584 0 4.56761e+10 1600 0 4.56788e+10 -8.8847e+06 5.09731e+06 -8.8847e+06 4.56587e+10 6.93062e+06 5.09731e+06 6.93062e+06 4.56659e+10 8.0177e-15 1.08148e+07 + + +4.26171e+06 2.4605e+06 0 -0.000689281 -0.112987 0 4.56467e+10 1600 -9.13492e-17 4.56524e+10 -9.73295e+06 7.29337e+06 -9.73295e+06 4.5644e+10 4.68681e+06 7.29337e+06 4.68681e+06 4.56313e+10 8.39786e-15 -1.85507e+07 +4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56771e+10 -2.06175e+07 -1.73077e+07 -2.06175e+07 4.56545e+10 -1.58625e+07 -1.73077e+07 -1.58625e+07 4.56938e+10 1.84777e-14 4.8163e+06 + + +5.51745e+06 3.1855e+06 0 -0.00516999 0.00895468 0 1.11727e+07 1600 1 -3.23725e+07 1.88972e+08 6.55746e+07 1.88972e+08 -9.6888e+08 5.55838e+07 6.55746e+07 5.55838e+07 2.56716e+08 3.36731e-16 1.11727e+07 +4.0952e+06 4.88047e+06 0 -0.063111 0.0529564 0 -1.38844e+07 1600 1 1.23456e+07 3.39785e+09 -2.70267e+08 3.39785e+09 1.07996e+07 2.20761e+08 -2.70267e+08 2.20761e+08 2.83632e+09 1.89255e-15 -1.38844e+07 + + +5.51745e+06 3.1855e+06 0 -0.00516999 0.00895468 0 1.11727e+07 1600 1 -3.23725e+07 1.88972e+08 6.55746e+07 1.88972e+08 -9.6888e+08 5.55838e+07 6.55746e+07 5.55838e+07 2.56716e+08 3.36731e-16 1.11727e+07 +5.18471e+06 2.99339e+06 2.17901e+06 -0.0171684 0.0586576 -0.0397298 -1.83017e+06 1600 1 3.28783e+09 3.56038e+09 2.14447e+09 3.56038e+09 2.49725e+07 -5.45396e+08 2.14447e+09 -5.45396e+08 2.94865e+09 2.27984e-15 -1.83017e+06 + + +3.16316e+06 3.7697e+06 0 0.0866503 -0.288584 0 4.56761e+10 1600 0 4.56788e+10 -8.8847e+06 5.09731e+06 -8.8847e+06 4.56587e+10 6.93062e+06 5.09731e+06 6.93062e+06 4.56659e+10 8.0177e-15 1.08148e+07 +4.0952e+06 4.88047e+06 0 -0.063111 0.0529564 0 -1.38844e+07 1600 1 1.23456e+07 3.39785e+09 -2.70267e+08 3.39785e+09 1.07996e+07 2.20761e+08 -2.70267e+08 2.20761e+08 2.83632e+09 1.89255e-15 -1.38844e+07 + + +3.16316e+06 3.7697e+06 0 0.0866503 -0.288584 0 4.56761e+10 1600 0 4.56788e+10 -8.8847e+06 5.09731e+06 -8.8847e+06 4.56587e+10 6.93062e+06 5.09731e+06 6.93062e+06 4.56659e+10 8.0177e-15 1.08148e+07 +2.9724e+06 3.54236e+06 1.68308e+06 -0.120017 -0.203298 -0.283249 4.56924e+10 1600 0 4.56611e+10 -7.03291e+06 -1.16048e+07 -7.03291e+06 4.56839e+10 -2.49822e+07 -1.16048e+07 -2.49822e+07 4.57141e+10 1.94695e-14 2.71349e+07 + + +4.0952e+06 4.88047e+06 0 -0.063111 0.0529564 0 -1.38844e+07 1600 1 1.23456e+07 3.39785e+09 -2.70267e+08 3.39785e+09 1.07996e+07 2.20761e+08 -2.70267e+08 2.20761e+08 2.83632e+09 1.89255e-15 -1.38844e+07 +3.84823e+06 4.58614e+06 2.17901e+06 -0.0834525 0.113859 -0.0922569 -1.96489e+07 1600 1 3.3168e+08 5.08411e+09 1.54245e+09 5.08411e+09 4.54904e+09 1.87102e+09 1.54245e+09 1.87102e+09 4.00823e+09 3.04113e-15 -1.96489e+07 + + +4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56771e+10 -2.06175e+07 -1.73077e+07 -2.06175e+07 4.56545e+10 -1.58625e+07 -1.73077e+07 -1.58625e+07 4.56938e+10 1.84777e-14 4.8163e+06 +5.18471e+06 2.99339e+06 2.17901e+06 -0.0171684 0.0586576 -0.0397298 -1.83017e+06 1600 1 3.28783e+09 3.56038e+09 2.14447e+09 3.56038e+09 2.49725e+07 -5.45396e+08 2.14447e+09 -5.45396e+08 2.94865e+09 2.27984e-15 -1.83017e+06 + + +4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56771e+10 -2.06175e+07 -1.73077e+07 -2.06175e+07 4.56545e+10 -1.58625e+07 -1.73077e+07 -1.58625e+07 4.56938e+10 1.84777e-14 4.8163e+06 +2.9724e+06 3.54236e+06 1.68308e+06 -0.120017 -0.203298 -0.283249 4.56924e+10 1600 0 4.56611e+10 -7.03291e+06 -1.16048e+07 -7.03291e+06 4.56839e+10 -2.49822e+07 -1.16048e+07 -2.49822e+07 4.57141e+10 1.94695e-14 2.71349e+07 + + +5.18471e+06 2.99339e+06 2.17901e+06 -0.0171684 0.0586576 -0.0397298 -1.83017e+06 1600 1 3.28783e+09 3.56038e+09 2.14447e+09 3.56038e+09 2.49725e+07 -5.45396e+08 2.14447e+09 -5.45396e+08 2.94865e+09 2.27984e-15 -1.83017e+06 +3.84823e+06 4.58614e+06 2.17901e+06 -0.0834525 0.113859 -0.0922569 -1.96489e+07 1600 1 3.3168e+08 5.08411e+09 1.54245e+09 5.08411e+09 4.54904e+09 1.87102e+09 1.54245e+09 1.87102e+09 4.00823e+09 3.04113e-15 -1.96489e+07 + + +2.9724e+06 3.54236e+06 1.68308e+06 -0.120017 -0.203298 -0.283249 4.56924e+10 1600 0 4.56611e+10 -7.03291e+06 -1.16048e+07 -7.03291e+06 4.56839e+10 -2.49822e+07 -1.16048e+07 -2.49822e+07 4.57141e+10 1.94695e-14 2.71349e+07 +3.84823e+06 4.58614e+06 2.17901e+06 -0.0834525 0.113859 -0.0922569 -1.96489e+07 1600 1 3.3168e+08 5.08411e+09 1.54245e+09 5.08411e+09 4.54904e+09 1.87102e+09 1.54245e+09 1.87102e+09 4.00823e+09 3.04113e-15 -1.96489e+07 + + +3.21212e+06 566384 1.18715e+06 0.185945 0.0327872 -0.518762 9.13123e+10 1600 0 9.12813e+10 -6.57581e+07 2.70157e+08 -6.57581e+07 9.12647e+10 2.31157e+08 2.70157e+08 2.31157e+08 9.05914e+10 2.67246e-13 -6.68253e+07 +4.55397e+06 802989 1.68308e+06 -0.148579 -0.0261985 -0.247956 4.56703e+10 1600 -5.52646e-17 4.56345e+10 -1.70458e+07 -1.33254e+06 -1.70458e+07 4.56638e+10 -43540.7 -1.33254e+06 -43540.7 4.56003e+10 1.804e-14 5.05121e+06 + + +3.21212e+06 566384 1.18715e+06 0.185945 0.0327872 -0.518762 9.13123e+10 1600 0 9.12813e+10 -6.57581e+07 2.70157e+08 -6.57581e+07 9.12647e+10 2.31157e+08 2.70157e+08 2.31157e+08 9.05914e+10 2.67246e-13 -6.68253e+07 +2.82469e+06 1.63084e+06 1.18715e+06 0.19034 0.0376886 -0.504666 9.13093e+10 1600 0 9.12035e+10 2.02277e+07 3.60459e+08 2.02277e+07 9.13395e+10 -8.36812e+06 3.60459e+08 -8.36812e+06 9.05849e+10 2.70251e-13 -6.9817e+07 + + +3.21212e+06 566384 1.18715e+06 0.185945 0.0327872 -0.518762 9.13123e+10 1600 0 9.12813e+10 -6.57581e+07 2.70157e+08 -6.57581e+07 9.12647e+10 2.31157e+08 2.70157e+08 2.31157e+08 9.05914e+10 2.67246e-13 -6.68253e+07 +2.61854e+06 461720 2.23112e+06 0 0 0 9.15676e+10 1600 0 9.18492e+10 4.86143e+07 -3.34847e+08 4.86143e+07 9.15779e+10 -5.90422e+07 -3.34847e+08 -5.90422e+07 9.21224e+10 2.19142e-13 1.88529e+08 + + +4.55397e+06 802989 1.68308e+06 -0.148579 -0.0261985 -0.247956 4.56703e+10 1600 -5.52646e-17 4.56345e+10 -1.70458e+07 -1.33254e+06 -1.70458e+07 4.56638e+10 -43540.7 -1.33254e+06 -43540.7 4.56003e+10 1.804e-14 5.05121e+06 +4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56315e+10 -1.35334e+07 1.17341e+06 -1.35334e+07 4.56694e+10 -4.32639e+06 1.17341e+06 -4.32639e+06 4.55995e+10 1.89028e-14 4.8163e+06 + + +4.55397e+06 802989 1.68308e+06 -0.148579 -0.0261985 -0.247956 4.56703e+10 1600 -5.52646e-17 4.56345e+10 -1.70458e+07 -1.33254e+06 -1.70458e+07 4.56638e+10 -43540.7 -1.33254e+06 -43540.7 4.56003e+10 1.804e-14 5.05121e+06 +3.71243e+06 654602 3.16316e+06 -0.708369 -0.124905 -0.603561 4.56608e+10 1600 -1.94027e-17 4.56527e+10 -1.19489e+08 2.62342e+07 -1.19489e+08 4.55834e+10 -1.0272e+08 2.62342e+07 -1.0272e+08 4.57566e+10 9.09935e-14 -4.52843e+06 + + +2.82469e+06 1.63084e+06 1.18715e+06 0.19034 0.0376886 -0.504666 9.13093e+10 1600 0 9.12035e+10 2.02277e+07 3.60459e+08 2.02277e+07 9.13395e+10 -8.36812e+06 3.60459e+08 -8.36812e+06 9.05849e+10 2.70251e-13 -6.9817e+07 +4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56315e+10 -1.35334e+07 1.17341e+06 -1.35334e+07 4.56694e+10 -4.32639e+06 1.17341e+06 -4.32639e+06 4.55995e+10 1.89028e-14 4.8163e+06 + + +2.82469e+06 1.63084e+06 1.18715e+06 0.19034 0.0376886 -0.504666 9.13093e+10 1600 0 9.12035e+10 2.02277e+07 3.60459e+08 2.02277e+07 9.13395e+10 -8.36812e+06 3.60459e+08 -8.36812e+06 9.05849e+10 2.70251e-13 -6.9817e+07 +2.30271e+06 1.32947e+06 2.23112e+06 -0.045879 0.0794648 0 9.15743e+10 1600 0 9.17898e+10 1.34251e+08 -2.95352e+08 1.34251e+08 9.16379e+10 -1.71012e+08 -2.95352e+08 -1.71012e+08 9.21342e+10 2.23142e-13 1.95234e+08 + + +4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56315e+10 -1.35334e+07 1.17341e+06 -1.35334e+07 4.56694e+10 -4.32639e+06 1.17341e+06 -4.32639e+06 4.55995e+10 1.89028e-14 4.8163e+06 +3.26466e+06 1.88485e+06 3.16316e+06 -0.659004 -0.297874 -0.603859 4.56587e+10 1600 -1.34824e-18 4.55015e+10 5.59097e+07 -4.75938e+07 5.59097e+07 4.57279e+10 9.34442e+07 -4.75938e+07 9.34442e+07 4.57547e+10 9.14777e-14 -6.62725e+06 + + +2.61854e+06 461720 2.23112e+06 0 0 0 9.15676e+10 1600 0 9.18492e+10 4.86143e+07 -3.34847e+08 4.86143e+07 9.15779e+10 -5.90422e+07 -3.34847e+08 -5.90422e+07 9.21224e+10 2.19142e-13 1.88529e+08 +3.71243e+06 654602 3.16316e+06 -0.708369 -0.124905 -0.603561 4.56608e+10 1600 -1.94027e-17 4.56527e+10 -1.19489e+08 2.62342e+07 -1.19489e+08 4.55834e+10 -1.0272e+08 2.62342e+07 -1.0272e+08 4.57566e+10 9.09935e-14 -4.52843e+06 + + +2.61854e+06 461720 2.23112e+06 0 0 0 9.15676e+10 1600 0 9.18492e+10 4.86143e+07 -3.34847e+08 4.86143e+07 9.15779e+10 -5.90422e+07 -3.34847e+08 -5.90422e+07 9.21224e+10 2.19142e-13 1.88529e+08 +2.30271e+06 1.32947e+06 2.23112e+06 -0.045879 0.0794648 0 9.15743e+10 1600 0 9.17898e+10 1.34251e+08 -2.95352e+08 1.34251e+08 9.16379e+10 -1.71012e+08 -2.95352e+08 -1.71012e+08 9.21342e+10 2.23142e-13 1.95234e+08 + + +3.71243e+06 654602 3.16316e+06 -0.708369 -0.124905 -0.603561 4.56608e+10 1600 -1.94027e-17 4.56527e+10 -1.19489e+08 2.62342e+07 -1.19489e+08 4.55834e+10 -1.0272e+08 2.62342e+07 -1.0272e+08 4.57566e+10 9.09935e-14 -4.52843e+06 +3.26466e+06 1.88485e+06 3.16316e+06 -0.659004 -0.297874 -0.603859 4.56587e+10 1600 -1.34824e-18 4.55015e+10 5.59097e+07 -4.75938e+07 5.59097e+07 4.57279e+10 9.34442e+07 -4.75938e+07 9.34442e+07 4.57547e+10 9.14777e-14 -6.62725e+06 + + +2.30271e+06 1.32947e+06 2.23112e+06 -0.045879 0.0794648 0 9.15743e+10 1600 0 9.17898e+10 1.34251e+08 -2.95352e+08 1.34251e+08 9.16379e+10 -1.71012e+08 -2.95352e+08 -1.71012e+08 9.21342e+10 2.23142e-13 1.95234e+08 +3.26466e+06 1.88485e+06 3.16316e+06 -0.659004 -0.297874 -0.603859 4.56587e+10 1600 -1.34824e-18 4.55015e+10 5.59097e+07 -4.75938e+07 5.59097e+07 4.57279e+10 9.34442e+07 -4.75938e+07 9.34442e+07 4.57547e+10 9.14777e-14 -6.62725e+06 + + +4.55397e+06 802989 1.68308e+06 -0.148579 -0.0261985 -0.247956 4.56703e+10 1600 -5.52646e-17 4.56737e+10 -1.01281e+07 -2.2518e+07 -1.01281e+07 4.5665e+10 -3.77906e+06 -2.2518e+07 -3.77906e+06 4.55793e+10 2.89328e-14 5.05121e+06 +5.89583e+06 1.03959e+06 2.17901e+06 0.0170632 0.00300871 -0.0476041 -3.4701e+06 1600 1 4.58898e+09 1.02622e+09 2.57071e+09 1.02622e+09 -1.18106e+09 3.64013e+08 2.57071e+09 3.64013e+08 -5.26011e+08 2.10803e-15 -3.4701e+06 + + +4.55397e+06 802989 1.68308e+06 -0.148579 -0.0261985 -0.247956 4.56703e+10 1600 -5.52646e-17 4.56737e+10 -1.01281e+07 -2.2518e+07 -1.01281e+07 4.5665e+10 -3.77906e+06 -2.2518e+07 -3.77906e+06 4.55793e+10 2.89328e-14 5.05121e+06 +4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.5662e+10 2.89722e+06 -1.71234e+07 2.89722e+06 4.56783e+10 -1.53802e+07 -1.71234e+07 -1.53802e+07 4.55787e+10 2.91045e-14 4.8163e+06 + + +4.55397e+06 802989 1.68308e+06 -0.148579 -0.0261985 -0.247956 4.56703e+10 1600 -5.52646e-17 4.56737e+10 -1.01281e+07 -2.2518e+07 -1.01281e+07 4.5665e+10 -3.77906e+06 -2.2518e+07 -3.77906e+06 4.55793e+10 2.89328e-14 5.05121e+06 +3.71243e+06 654602 3.16316e+06 -0.708369 -0.124905 -0.603561 4.56608e+10 1600 -1.94027e-17 4.5615e+10 -1.26129e+08 -5.85141e+06 -1.26129e+08 4.55822e+10 -1.08378e+08 -5.85141e+06 -1.08378e+08 4.57293e+10 9.17168e-14 -4.52843e+06 + + +5.89583e+06 1.03959e+06 2.17901e+06 0.0170632 0.00300871 -0.0476041 -3.4701e+06 1600 1 4.58898e+09 1.02622e+09 2.57071e+09 1.02622e+09 -1.18106e+09 3.64013e+08 2.57071e+09 3.64013e+08 -5.26011e+08 2.10803e-15 -3.4701e+06 +5.18471e+06 2.99339e+06 2.17901e+06 -0.0171684 0.0586576 -0.0397298 -1.83017e+06 1600 1 2.74028e+09 3.63106e+09 3.43815e+09 3.63106e+09 -3.28105e+08 -2.93871e+08 3.43815e+09 -2.93871e+08 -5.94477e+08 2.6705e-15 -1.83017e+06 + + +5.89583e+06 1.03959e+06 2.17901e+06 0.0170632 0.00300871 -0.0476041 -3.4701e+06 1600 1 4.58898e+09 1.02622e+09 2.57071e+09 1.02622e+09 -1.18106e+09 3.64013e+08 2.57071e+09 3.64013e+08 -5.26011e+08 2.10803e-15 -3.4701e+06 +4.80632e+06 847485 4.0952e+06 0 0 0 1.81123e+07 1600 1 -1.96764e+10 -2.92287e+09 -1.55519e+10 -2.92287e+09 -4.19079e+09 -2.74221e+09 -1.55519e+10 -2.74221e+09 -1.46328e+10 8.94851e-15 1.81123e+07 + + +4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.5662e+10 2.89722e+06 -1.71234e+07 2.89722e+06 4.56783e+10 -1.53802e+07 -1.71234e+07 -1.53802e+07 4.55787e+10 2.91045e-14 4.8163e+06 +5.18471e+06 2.99339e+06 2.17901e+06 -0.0171684 0.0586576 -0.0397298 -1.83017e+06 1600 1 2.74028e+09 3.63106e+09 3.43815e+09 3.63106e+09 -3.28105e+08 -2.93871e+08 3.43815e+09 -2.93871e+08 -5.94477e+08 2.6705e-15 -1.83017e+06 + + +4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.5662e+10 2.89722e+06 -1.71234e+07 2.89722e+06 4.56783e+10 -1.53802e+07 -1.71234e+07 -1.53802e+07 4.55787e+10 2.91045e-14 4.8163e+06 +3.26466e+06 1.88485e+06 3.16316e+06 -0.659004 -0.297874 -0.603859 4.56587e+10 1600 -1.34824e-18 4.54752e+10 3.89659e+07 -7.38383e+07 3.89659e+07 4.57171e+10 7.66056e+07 -7.38383e+07 7.66056e+07 4.57286e+10 9.127e-14 -6.62725e+06 + + +5.18471e+06 2.99339e+06 2.17901e+06 -0.0171684 0.0586576 -0.0397298 -1.83017e+06 1600 1 2.74028e+09 3.63106e+09 3.43815e+09 3.63106e+09 -3.28105e+08 -2.93871e+08 3.43815e+09 -2.93871e+08 -5.94477e+08 2.6705e-15 -1.83017e+06 +4.22661e+06 2.44023e+06 4.0952e+06 -0.0644111 0.111563 3.93309e-17 1.73775e+07 1600 1 -1.68018e+10 -5.75744e+09 -1.36291e+10 -5.75744e+09 -9.28812e+09 -8.51327e+09 -1.36291e+10 -8.51327e+09 -1.41625e+10 8.74508e-15 1.73775e+07 + + +3.71243e+06 654602 3.16316e+06 -0.708369 -0.124905 -0.603561 4.56608e+10 1600 -1.94027e-17 4.5615e+10 -1.26129e+08 -5.85141e+06 -1.26129e+08 4.55822e+10 -1.08378e+08 -5.85141e+06 -1.08378e+08 4.57293e+10 9.17168e-14 -4.52843e+06 +4.80632e+06 847485 4.0952e+06 0 0 0 1.81123e+07 1600 1 -1.96764e+10 -2.92287e+09 -1.55519e+10 -2.92287e+09 -4.19079e+09 -2.74221e+09 -1.55519e+10 -2.74221e+09 -1.46328e+10 8.94851e-15 1.81123e+07 + + +3.71243e+06 654602 3.16316e+06 -0.708369 -0.124905 -0.603561 4.56608e+10 1600 -1.94027e-17 4.5615e+10 -1.26129e+08 -5.85141e+06 -1.26129e+08 4.55822e+10 -1.08378e+08 -5.85141e+06 -1.08378e+08 4.57293e+10 9.17168e-14 -4.52843e+06 +3.26466e+06 1.88485e+06 3.16316e+06 -0.659004 -0.297874 -0.603859 4.56587e+10 1600 -1.34824e-18 4.54752e+10 3.89659e+07 -7.38383e+07 3.89659e+07 4.57171e+10 7.66056e+07 -7.38383e+07 7.66056e+07 4.57286e+10 9.127e-14 -6.62725e+06 + + +4.80632e+06 847485 4.0952e+06 0 0 0 1.81123e+07 1600 1 -1.96764e+10 -2.92287e+09 -1.55519e+10 -2.92287e+09 -4.19079e+09 -2.74221e+09 -1.55519e+10 -2.74221e+09 -1.46328e+10 8.94851e-15 1.81123e+07 +4.22661e+06 2.44023e+06 4.0952e+06 -0.0644111 0.111563 3.93309e-17 1.73775e+07 1600 1 -1.68018e+10 -5.75744e+09 -1.36291e+10 -5.75744e+09 -9.28812e+09 -8.51327e+09 -1.36291e+10 -8.51327e+09 -1.41625e+10 8.74508e-15 1.73775e+07 + + +3.26466e+06 1.88485e+06 3.16316e+06 -0.659004 -0.297874 -0.603859 4.56587e+10 1600 -1.34824e-18 4.54752e+10 3.89659e+07 -7.38383e+07 3.89659e+07 4.57171e+10 7.66056e+07 -7.38383e+07 7.66056e+07 4.57286e+10 9.127e-14 -6.62725e+06 +4.22661e+06 2.44023e+06 4.0952e+06 -0.0644111 0.111563 3.93309e-17 1.73775e+07 1600 1 -1.68018e+10 -5.75744e+09 -1.36291e+10 -5.75744e+09 -9.28812e+09 -8.51327e+09 -1.36291e+10 -8.51327e+09 -1.41625e+10 8.74508e-15 1.73775e+07 + + +2.82469e+06 1.63084e+06 1.18715e+06 0.19034 0.0376886 -0.504666 9.13093e+10 1600 0 9.13174e+10 -4.73794e+07 1.8135e+08 -4.73794e+07 9.12319e+10 3.01858e+08 1.8135e+08 3.01858e+08 9.05849e+10 2.67694e-13 -6.9817e+07 +4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56578e+10 -3.05214e+07 1.47808e+06 -3.05214e+07 4.56495e+10 -4.8541e+06 1.47808e+06 -4.8541e+06 4.55995e+10 2.20884e-14 4.8163e+06 + + +2.82469e+06 1.63084e+06 1.18715e+06 0.19034 0.0376886 -0.504666 9.13093e+10 1600 0 9.13174e+10 -4.73794e+07 1.8135e+08 -4.73794e+07 9.12319e+10 3.01858e+08 1.8135e+08 3.01858e+08 9.05849e+10 2.67694e-13 -6.9817e+07 +2.09656e+06 2.49859e+06 1.18715e+06 0.151969 0.147411 -0.578639 9.13366e+10 1600 0 9.11902e+10 -2.31881e+07 3.58683e+08 -2.31881e+07 9.13885e+10 1.07386e+08 3.58683e+08 1.07386e+08 9.05995e+10 2.78025e-13 -4.25455e+07 + + +2.82469e+06 1.63084e+06 1.18715e+06 0.19034 0.0376886 -0.504666 9.13093e+10 1600 0 9.13174e+10 -4.73794e+07 1.8135e+08 -4.73794e+07 9.12319e+10 3.01858e+08 1.8135e+08 3.01858e+08 9.05849e+10 2.67694e-13 -6.9817e+07 +2.30271e+06 1.32947e+06 2.23112e+06 -0.045879 0.0794648 0 9.15743e+10 1600 0 9.17795e+10 1.47943e+08 -2.95352e+08 1.47943e+08 9.16213e+10 -1.71012e+08 -2.95352e+08 -1.71012e+08 9.21342e+10 2.27683e-13 1.95234e+08 + + +4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56578e+10 -3.05214e+07 1.47808e+06 -3.05214e+07 4.56495e+10 -4.8541e+06 1.47808e+06 -4.8541e+06 4.55995e+10 2.20884e-14 4.8163e+06 +2.9724e+06 3.54236e+06 1.68308e+06 -0.120017 -0.203298 -0.283249 4.56924e+10 1600 0 4.56808e+10 -7.47783e+06 -5.79055e+06 -7.47783e+06 4.56548e+10 183550 -5.79055e+06 183550 4.56234e+10 1.51229e-14 2.71349e+07 + + +4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56578e+10 -3.05214e+07 1.47808e+06 -3.05214e+07 4.56495e+10 -4.8541e+06 1.47808e+06 -4.8541e+06 4.55995e+10 2.20884e-14 4.8163e+06 +3.26466e+06 1.88485e+06 3.16316e+06 -0.659004 -0.297874 -0.603859 4.56587e+10 1600 -1.34824e-18 4.57179e+10 -7.08523e+07 5.76269e+07 -7.08523e+07 4.55179e+10 -8.88035e+07 5.76269e+07 -8.88035e+07 4.57547e+10 9.00879e-14 -6.62725e+06 + + +2.09656e+06 2.49859e+06 1.18715e+06 0.151969 0.147411 -0.578639 9.13366e+10 1600 0 9.11902e+10 -2.31881e+07 3.58683e+08 -2.31881e+07 9.13885e+10 1.07386e+08 3.58683e+08 1.07386e+08 9.05995e+10 2.78025e-13 -4.25455e+07 +2.9724e+06 3.54236e+06 1.68308e+06 -0.120017 -0.203298 -0.283249 4.56924e+10 1600 0 4.56808e+10 -7.47783e+06 -5.79055e+06 -7.47783e+06 4.56548e+10 183550 -5.79055e+06 183550 4.56234e+10 1.51229e-14 2.71349e+07 + + +2.09656e+06 2.49859e+06 1.18715e+06 0.151969 0.147411 -0.578639 9.13366e+10 1600 0 9.11902e+10 -2.31881e+07 3.58683e+08 -2.31881e+07 9.13885e+10 1.07386e+08 3.58683e+08 1.07386e+08 9.05995e+10 2.78025e-13 -4.25455e+07 +1.70913e+06 2.03687e+06 2.23112e+06 -0.11234 0.0942646 0 9.15041e+10 1600 0 9.16788e+10 1.36191e+08 -2.31722e+08 1.36191e+08 9.16625e+10 -2.50142e+08 -2.31722e+08 -2.50142e+08 9.2071e+10 2.16969e-13 1.24956e+08 + + +2.9724e+06 3.54236e+06 1.68308e+06 -0.120017 -0.203298 -0.283249 4.56924e+10 1600 0 4.56808e+10 -7.47783e+06 -5.79055e+06 -7.47783e+06 4.56548e+10 183550 -5.79055e+06 183550 4.56234e+10 1.51229e-14 2.71349e+07 +2.42312e+06 2.88776e+06 3.16316e+06 -0.690435 -0.372828 -0.612044 4.56604e+10 1600 0 4.54595e+10 -2.64945e+07 -6.95638e+07 -2.64945e+07 4.57508e+10 6.4528e+07 -6.95638e+07 6.4528e+07 4.575e+10 9.73365e-14 -4.90659e+06 + + +2.30271e+06 1.32947e+06 2.23112e+06 -0.045879 0.0794648 0 9.15743e+10 1600 0 9.17795e+10 1.47943e+08 -2.95352e+08 1.47943e+08 9.16213e+10 -1.71012e+08 -2.95352e+08 -1.71012e+08 9.21342e+10 2.27683e-13 1.95234e+08 +3.26466e+06 1.88485e+06 3.16316e+06 -0.659004 -0.297874 -0.603859 4.56587e+10 1600 -1.34824e-18 4.57179e+10 -7.08523e+07 5.76269e+07 -7.08523e+07 4.55179e+10 -8.88035e+07 5.76269e+07 -8.88035e+07 4.57547e+10 9.00879e-14 -6.62725e+06 + + +2.30271e+06 1.32947e+06 2.23112e+06 -0.045879 0.0794648 0 9.15743e+10 1600 0 9.17795e+10 1.47943e+08 -2.95352e+08 1.47943e+08 9.16213e+10 -1.71012e+08 -2.95352e+08 -1.71012e+08 9.21342e+10 2.27683e-13 1.95234e+08 +1.70913e+06 2.03687e+06 2.23112e+06 -0.11234 0.0942646 0 9.15041e+10 1600 0 9.16788e+10 1.36191e+08 -2.31722e+08 1.36191e+08 9.16625e+10 -2.50142e+08 -2.31722e+08 -2.50142e+08 9.2071e+10 2.16969e-13 1.24956e+08 + + +3.26466e+06 1.88485e+06 3.16316e+06 -0.659004 -0.297874 -0.603859 4.56587e+10 1600 -1.34824e-18 4.57179e+10 -7.08523e+07 5.76269e+07 -7.08523e+07 4.55179e+10 -8.88035e+07 5.76269e+07 -8.88035e+07 4.57547e+10 9.00879e-14 -6.62725e+06 +2.42312e+06 2.88776e+06 3.16316e+06 -0.690435 -0.372828 -0.612044 4.56604e+10 1600 0 4.54595e+10 -2.64945e+07 -6.95638e+07 -2.64945e+07 4.57508e+10 6.4528e+07 -6.95638e+07 6.4528e+07 4.575e+10 9.73365e-14 -4.90659e+06 + + +1.70913e+06 2.03687e+06 2.23112e+06 -0.11234 0.0942646 0 9.15041e+10 1600 0 9.16788e+10 1.36191e+08 -2.31722e+08 1.36191e+08 9.16625e+10 -2.50142e+08 -2.31722e+08 -2.50142e+08 9.2071e+10 2.16969e-13 1.24956e+08 +2.42312e+06 2.88776e+06 3.16316e+06 -0.690435 -0.372828 -0.612044 4.56604e+10 1600 0 4.54595e+10 -2.64945e+07 -6.95638e+07 -2.64945e+07 4.57508e+10 6.4528e+07 -6.95638e+07 6.4528e+07 4.575e+10 9.73365e-14 -4.90659e+06 + + +4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56882e+10 -1.40908e+07 -1.68188e+07 -1.40908e+07 4.56583e+10 -1.59078e+07 -1.68188e+07 -1.59078e+07 4.55787e+10 3.13788e-14 4.8163e+06 +5.18471e+06 2.99339e+06 2.17901e+06 -0.0171684 0.0586576 -0.0397298 -1.83017e+06 1600 1 2.85262e+09 3.43257e+09 3.39323e+09 3.43257e+09 2.24543e+07 -2.1608e+08 3.39323e+09 -2.1608e+08 -5.94477e+08 2.58467e-15 -1.83017e+06 + + +4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56882e+10 -1.40908e+07 -1.68188e+07 -1.40908e+07 4.56583e+10 -1.59078e+07 -1.68188e+07 -1.59078e+07 4.55787e+10 3.13788e-14 4.8163e+06 +2.9724e+06 3.54236e+06 1.68308e+06 -0.120017 -0.203298 -0.283249 4.56924e+10 1600 0 4.56741e+10 1.9163e+06 -2.63568e+07 1.9163e+06 4.56867e+10 -1.44764e+07 -2.63568e+07 -1.44764e+07 4.56023e+10 2.72956e-14 2.71349e+07 + + +4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56882e+10 -1.40908e+07 -1.68188e+07 -1.40908e+07 4.56583e+10 -1.59078e+07 -1.68188e+07 -1.59078e+07 4.55787e+10 3.13788e-14 4.8163e+06 +3.26466e+06 1.88485e+06 3.16316e+06 -0.659004 -0.297874 -0.603859 4.56587e+10 1600 -1.34824e-18 4.56915e+10 -8.7796e+07 3.13824e+07 -8.7796e+07 4.55071e+10 -1.05642e+08 3.13824e+07 -1.05642e+08 4.57286e+10 9.20921e-14 -6.62725e+06 + + +5.18471e+06 2.99339e+06 2.17901e+06 -0.0171684 0.0586576 -0.0397298 -1.83017e+06 1600 1 2.85262e+09 3.43257e+09 3.39323e+09 3.43257e+09 2.24543e+07 -2.1608e+08 3.39323e+09 -2.1608e+08 -5.94477e+08 2.58467e-15 -1.83017e+06 +3.84823e+06 4.58614e+06 2.17901e+06 -0.0834525 0.113859 -0.0922569 -1.96489e+07 1600 1 4.76835e+08 5.04687e+09 1.53814e+09 5.04687e+09 4.25411e+09 2.76449e+09 1.53814e+09 2.76449e+09 1.39318e+09 3.13692e-15 -1.96489e+07 + + +5.18471e+06 2.99339e+06 2.17901e+06 -0.0171684 0.0586576 -0.0397298 -1.83017e+06 1600 1 2.85262e+09 3.43257e+09 3.39323e+09 3.43257e+09 2.24543e+07 -2.1608e+08 3.39323e+09 -2.1608e+08 -5.94477e+08 2.58467e-15 -1.83017e+06 +4.22661e+06 2.44023e+06 4.0952e+06 -0.0644111 0.111563 3.93309e-17 1.73775e+07 1600 1 -1.6204e+10 -6.51291e+09 -1.36291e+10 -6.51291e+09 -8.46455e+09 -8.51325e+09 -1.36291e+10 -8.51325e+09 -1.41625e+10 8.89851e-15 1.73775e+07 + + +2.9724e+06 3.54236e+06 1.68308e+06 -0.120017 -0.203298 -0.283249 4.56924e+10 1600 0 4.56741e+10 1.9163e+06 -2.63568e+07 1.9163e+06 4.56867e+10 -1.44764e+07 -2.63568e+07 -1.44764e+07 4.56023e+10 2.72956e-14 2.71349e+07 +3.84823e+06 4.58614e+06 2.17901e+06 -0.0834525 0.113859 -0.0922569 -1.96489e+07 1600 1 4.76835e+08 5.04687e+09 1.53814e+09 5.04687e+09 4.25411e+09 2.76449e+09 1.53814e+09 2.76449e+09 1.39318e+09 3.13692e-15 -1.96489e+07 + + +2.9724e+06 3.54236e+06 1.68308e+06 -0.120017 -0.203298 -0.283249 4.56924e+10 1600 0 4.56741e+10 1.9163e+06 -2.63568e+07 1.9163e+06 4.56867e+10 -1.44764e+07 -2.63568e+07 -1.44764e+07 4.56023e+10 2.72956e-14 2.71349e+07 +2.42312e+06 2.88776e+06 3.16316e+06 -0.690435 -0.372828 -0.612044 4.56604e+10 1600 0 4.54514e+10 -4.2775e+07 -8.43876e+07 -4.2775e+07 4.57234e+10 3.81818e+07 -8.43876e+07 3.81818e+07 4.5725e+10 9.38237e-14 -4.90659e+06 + + +3.84823e+06 4.58614e+06 2.17901e+06 -0.0834525 0.113859 -0.0922569 -1.96489e+07 1600 1 4.76835e+08 5.04687e+09 1.53814e+09 5.04687e+09 4.25411e+09 2.76449e+09 1.53814e+09 2.76449e+09 1.39318e+09 3.13692e-15 -1.96489e+07 +3.13711e+06 3.73866e+06 4.0952e+06 -0.291372 0.244491 2.21699e-16 -1.06766e+07 1600 1 -1.77671e+10 -3.19649e+09 -1.66819e+09 -3.19649e+09 -1.75129e+10 -1.52152e+10 -1.66819e+09 -1.52152e+10 -1.68353e+10 7.82202e-15 -1.06766e+07 + + +3.26466e+06 1.88485e+06 3.16316e+06 -0.659004 -0.297874 -0.603859 4.56587e+10 1600 -1.34824e-18 4.56915e+10 -8.7796e+07 3.13824e+07 -8.7796e+07 4.55071e+10 -1.05642e+08 3.13824e+07 -1.05642e+08 4.57286e+10 9.20921e-14 -6.62725e+06 +4.22661e+06 2.44023e+06 4.0952e+06 -0.0644111 0.111563 3.93309e-17 1.73775e+07 1600 1 -1.6204e+10 -6.51291e+09 -1.36291e+10 -6.51291e+09 -8.46455e+09 -8.51325e+09 -1.36291e+10 -8.51325e+09 -1.41625e+10 8.89851e-15 1.73775e+07 + + +3.26466e+06 1.88485e+06 3.16316e+06 -0.659004 -0.297874 -0.603859 4.56587e+10 1600 -1.34824e-18 4.56915e+10 -8.7796e+07 3.13824e+07 -8.7796e+07 4.55071e+10 -1.05642e+08 3.13824e+07 -1.05642e+08 4.57286e+10 9.20921e-14 -6.62725e+06 +2.42312e+06 2.88776e+06 3.16316e+06 -0.690435 -0.372828 -0.612044 4.56604e+10 1600 0 4.54514e+10 -4.2775e+07 -8.43876e+07 -4.2775e+07 4.57234e+10 3.81818e+07 -8.43876e+07 3.81818e+07 4.5725e+10 9.38237e-14 -4.90659e+06 + + +4.22661e+06 2.44023e+06 4.0952e+06 -0.0644111 0.111563 3.93309e-17 1.73775e+07 1600 1 -1.6204e+10 -6.51291e+09 -1.36291e+10 -6.51291e+09 -8.46455e+09 -8.51325e+09 -1.36291e+10 -8.51325e+09 -1.41625e+10 8.89851e-15 1.73775e+07 +3.13711e+06 3.73866e+06 4.0952e+06 -0.291372 0.244491 2.21699e-16 -1.06766e+07 1600 1 -1.77671e+10 -3.19649e+09 -1.66819e+09 -3.19649e+09 -1.75129e+10 -1.52152e+10 -1.66819e+09 -1.52152e+10 -1.68353e+10 7.82202e-15 -1.06766e+07 + + +2.42312e+06 2.88776e+06 3.16316e+06 -0.690435 -0.372828 -0.612044 4.56604e+10 1600 0 4.54514e+10 -4.2775e+07 -8.43876e+07 -4.2775e+07 4.57234e+10 3.81818e+07 -8.43876e+07 3.81818e+07 4.5725e+10 9.38237e-14 -4.90659e+06 +3.13711e+06 3.73866e+06 4.0952e+06 -0.291372 0.244491 2.21699e-16 -1.06766e+07 1600 1 -1.77671e+10 -3.19649e+09 -1.66819e+09 -3.19649e+09 -1.75129e+10 -1.52152e+10 -1.66819e+09 -1.52152e+10 -1.68353e+10 7.82202e-15 -1.06766e+07 + + diff --git a/tests/chunk_lithostatic_pressure_2d_initial_topography/statistics b/tests/chunk_lithostatic_pressure_2d_initial_topography/statistics new file mode 100644 index 00000000000..7c4ebdcbbd0 --- /dev/null +++ b/tests/chunk_lithostatic_pressure_2d_initial_topography/statistics @@ -0,0 +1,19 @@ +# 1: Time step number +# 2: Time (years) +# 3: Time step size (years) +# 4: Number of mesh cells +# 5: Number of Stokes degrees of freedom +# 6: Number of temperature degrees of freedom +# 7: Number of degrees of freedom for all compositions +# 8: Iterations for temperature solver +# 9: Iterations for composition solver 1 +# 10: Iterations for Stokes solver +# 11: Velocity iterations in Stokes preconditioner +# 12: Schur complement iterations in Stokes preconditioner +# 13: RMS velocity (m/year) +# 14: Max. velocity (m/year) +# 15: Minimal pressure (Pa) +# 16: Average pressure (Pa) +# 17: Maximal pressure (Pa) +# 18: Visualization file name +0 0.000000000000e+00 0.000000000000e+00 8 402 125 125 0 11 21 23 22 3.09329044e-01 1.86122381e+00 -1.96488886e+07 3.29000931e+10 9.15743340e+10 output-chunk_lithostatic_pressure_2d/solution/solution-00000 From 6ea426fd5d2a6b3af0e9542ab0f677eaeea59595 Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Fri, 8 Dec 2023 17:59:21 +0100 Subject: [PATCH 12/22] Test function over time --- tests/airy_isostasy_initial_topo.cc | 76 ++++++++++++++++++++++++++++ tests/airy_isostasy_initial_topo.prm | 5 ++ 2 files changed, 81 insertions(+) create mode 100644 tests/airy_isostasy_initial_topo.cc diff --git a/tests/airy_isostasy_initial_topo.cc b/tests/airy_isostasy_initial_topo.cc new file mode 100644 index 00000000000..9ad525be681 --- /dev/null +++ b/tests/airy_isostasy_initial_topo.cc @@ -0,0 +1,76 @@ +/* + Copyright (C) 2022 by the authors of the ASPECT code. + + This file is part of ASPECT. + + ASPECT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2, or (at your option) + any later version. + + ASPECT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ASPECT; see the file LICENSE. If not see + . +*/ + +#include +#include + + +namespace aspect +{ + namespace Postprocess + { + template + class PointIsInDomain : public Interface, public ::aspect::SimulatorAccess + { + public: + virtual + std::pair + execute (TableHandler &); + }; + } +} + + +namespace aspect +{ + namespace Postprocess + { + template + std::pair + PointIsInDomain::execute(TableHandler &) + { + // The airy isostasy box domain is 1x1 with max 0.05 topography. + const Point point_in_domain(0.1, 0.1); + const Point point_not_in_domain(0.1, 2.0); + const Point point_not_in_domain_2(2.0, 0.01); + + std::ostringstream screen_text; + screen_text << std::boolalpha << this->get_geometry_model().point_is_in_domain(point_in_domain) ; + screen_text << std::boolalpha << this->get_geometry_model().point_is_in_domain(point_not_in_domain); + screen_text << std::boolalpha << this->get_geometry_model().point_is_in_domain(point_not_in_domain_2); + + return std::pair ("Points lie in domain: ", + screen_text.str()); + } + } +} + + +// explicit instantiations +namespace aspect +{ + namespace Postprocess + { + ASPECT_REGISTER_POSTPROCESSOR(PointIsInDomain, + "point is in domain", + "A postprocessor that tests whether certain points " + "lie within the domain.") + } +} diff --git a/tests/airy_isostasy_initial_topo.prm b/tests/airy_isostasy_initial_topo.prm index f0a9420d49d..867ab8cb2e1 100644 --- a/tests/airy_isostasy_initial_topo.prm +++ b/tests/airy_isostasy_initial_topo.prm @@ -30,3 +30,8 @@ subsection Boundary traction model set Function expression = 0;10 end end + +# Also test the point_is_in_domain function of the geometry model. +subsection Postprocess + set List of postprocessors = velocity statistics, pressure statistics, topography, point is in domain +end From 4e976a5b227ee893e0a90963b660e2a2047fe8ff Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Fri, 8 Dec 2023 17:59:38 +0100 Subject: [PATCH 13/22] Update comment --- source/geometry_model/box.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/geometry_model/box.cc b/source/geometry_model/box.cc index 40045472929..c24efe31cce 100644 --- a/source/geometry_model/box.cc +++ b/source/geometry_model/box.cc @@ -354,7 +354,7 @@ namespace aspect Point max_point = extents+box_origin; // If mesh deformation is not enabled, but initial topography - // was applied to the mesh, include this topography in the + // was/will be applied to the mesh, include this topography in the // extent of the domain. if (!Plugins::plugin_type_matches>(this->get_initial_topography_model())) { From 8d96d2d4f18779ca215a0312058173dc14d0d566 Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Mon, 11 Dec 2023 17:44:51 +0100 Subject: [PATCH 14/22] Check all procs for active cell --- source/geometry_model/box.cc | 59 +++++++++++++++++++++------------- source/geometry_model/chunk.cc | 32 ++++++++++++------ 2 files changed, 59 insertions(+), 32 deletions(-) diff --git a/source/geometry_model/box.cc b/source/geometry_model/box.cc index c24efe31cce..6d17441884a 100644 --- a/source/geometry_model/box.cc +++ b/source/geometry_model/box.cc @@ -340,41 +340,54 @@ namespace aspect if (this->get_parameters().mesh_deformation_enabled && this->simulator_is_past_initialization()) { + bool cell_found = false; std::pair::active_cell_iterator, Point> it = GridTools::find_active_cell_around_point<> (this->get_mapping(), this->get_triangulation(), point); - if (it.first.state() == IteratorState::valid) + // If the cell the point is in is on this processor, we have found the right cell. + if (it.first.state() == IteratorState::valid && it.first->is_locally_owned()) + cell_found = true; + + // Compute how many processes found the cell. + const int n_procs_cell_found = Utilities::MPI::sum(cell_found ? 1 : 0, this->get_mpi_communicator()); + + // If at least one process found the cell, the point is in the domain. + if (n_procs_cell_found > 0) return true; else return false; } + // Without mesh deformation enabled, it is much cheaper to check whether the point lies in the domain. + else + { - // The maximal extents of the unperturbed box domain. - Point max_point = extents+box_origin; + // The maximal extents of the unperturbed box domain. + Point max_point = extents+box_origin; - // If mesh deformation is not enabled, but initial topography - // was/will be applied to the mesh, include this topography in the - // extent of the domain. - if (!Plugins::plugin_type_matches>(this->get_initial_topography_model())) - { - // Get the surface x (,y) point - Point surface_point; - for (unsigned int d=0; dvalue(surface_point); - max_point[dim-1] += topo; - } + // If mesh deformation is not enabled, but initial topography + // was/will be applied to the mesh, include this topography in the + // extent of the domain. + if (!Plugins::plugin_type_matches>(this->get_initial_topography_model())) + { + // Get the surface x (,y) point + Point surface_point; + for (unsigned int d=0; dvalue(surface_point); + max_point[dim-1] += topo; + } - // Check whether point lies within the min/max coordinates of the domain including initial topography. - for (unsigned int d = 0; d < dim; ++d) - if (point[d] > max_point[d]+std::numeric_limits::epsilon()*extents[d] || - point[d] < box_origin[d]-std::numeric_limits::epsilon()*extents[d]) - return false; + // Check whether point lies within the min/max coordinates of the domain including initial topography. + for (unsigned int d = 0; d < dim; ++d) + if (point[d] > max_point[d]+std::numeric_limits::epsilon()*extents[d] || + point[d] < box_origin[d]-std::numeric_limits::epsilon()*extents[d]) + return false; - return true; + return true; + } } template diff --git a/source/geometry_model/chunk.cc b/source/geometry_model/chunk.cc index 4fb70dec3eb..36049380d0c 100644 --- a/source/geometry_model/chunk.cc +++ b/source/geometry_model/chunk.cc @@ -707,23 +707,37 @@ namespace aspect if (this->get_parameters().mesh_deformation_enabled && this->simulator_is_past_initialization()) { + bool cell_found = false; std::pair::active_cell_iterator, - Point> it = - GridTools::find_active_cell_around_point<> (this->get_mapping(), this->get_triangulation(), point); + Point> + it = + GridTools::find_active_cell_around_point<>(this->get_mapping(), this->get_triangulation(), point); - if (it.first.state() == IteratorState::valid) + // If the cell the point is in is on this processor, we have found the right cell. + if (it.first.state() == IteratorState::valid && it.first->is_locally_owned()) + cell_found = true; + + // Compute how many processes found the cell. + const int n_procs_cell_found = Utilities::MPI::sum(cell_found ? 1 : 0, this->get_mpi_communicator()); + + // If at least one process found the cell, the point is in the domain. + if (n_procs_cell_found > 0) return true; else return false; } - const Point spherical_point = manifold->pull_back(point); + // Without mesh deformation enabled, it is much cheaper to check whether the point lies in the domain. + else + { + const Point spherical_point = manifold->pull_back(point); - for (unsigned int d = 0; d < dim; ++d) - if (spherical_point[d] > point2[d]+std::numeric_limits::epsilon()*std::abs(point2[d]) || - spherical_point[d] < point1[d]-std::numeric_limits::epsilon()*std::abs(point2[d])) - return false; + for (unsigned int d = 0; d < dim; ++d) + if (spherical_point[d] > point2[d]+std::numeric_limits::epsilon()*std::abs(point2[d]) || + spherical_point[d] < point1[d]-std::numeric_limits::epsilon()*std::abs(point2[d])) + return false; - return true; + return true; + } } From 7577387749e223503edc350c59475feccdde1783 Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Mon, 11 Dec 2023 22:03:18 +0100 Subject: [PATCH 15/22] Fix tests --- tests/airy_isostasy_initial_topo.cc | 1 + ...ostatic_pressure_2d_initial_topography.prm | 23 ++++++++++--------- 2 files changed, 13 insertions(+), 11 deletions(-) diff --git a/tests/airy_isostasy_initial_topo.cc b/tests/airy_isostasy_initial_topo.cc index 9ad525be681..3f194ab00e9 100644 --- a/tests/airy_isostasy_initial_topo.cc +++ b/tests/airy_isostasy_initial_topo.cc @@ -19,6 +19,7 @@ */ #include +#include #include diff --git a/tests/chunk_lithostatic_pressure_2d_initial_topography.prm b/tests/chunk_lithostatic_pressure_2d_initial_topography.prm index 7861c6a79ae..8f7243b24fb 100644 --- a/tests/chunk_lithostatic_pressure_2d_initial_topography.prm +++ b/tests/chunk_lithostatic_pressure_2d_initial_topography.prm @@ -8,7 +8,7 @@ include $ASPECT_SOURCE_DIR/tests/chunk_lithostatic_pressure_2d.prm -set Dimension = 3 +set Dimension = 2 set End time = 0 subsection Geometry model @@ -17,29 +17,30 @@ subsection Geometry model subsection Chunk set Chunk inner radius = 3471000 set Chunk outer radius = 6371000 - set Chunk minimum longitude = 10 + set Chunk minimum longitude = 0 set Chunk maximum longitude = 50 - set Chunk minimum latitude = 0 - set Chunk maximum latitude = 40 end subsection Initial topography model - set Model name = function + set Model name = ascii data - subsection Function - set Function expression = 5000 - set Maximum topography value = 5000 + subsection Ascii data model + set Data directory = $ASPECT_SOURCE_DIR/data/geometry-model/initial-topography-model/ascii-data/test/ + set Data file name = shell_2d_outer.0.txt end end end +subsection Boundary velocity model + set Tangential velocity boundary indicators = inner, outer, west +end subsection Boundary traction model set Prescribed traction boundary indicators = east: initial lithostatic pressure - # Set a radius between 6371 and 6376 km (surface radius - # without and with initial topography. + # Set a radius between 6371 and 6461 km (surface radius + # without and with 90 km initial topography at 0 degrees longitude. subsection Initial lithostatic pressure - set Representative point = 6374000, 50, 40 + set Representative point = 6394000, 0.1 end end From 8532588e153203f051085b755d6ea822405407a7 Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Thu, 14 Dec 2023 11:41:52 +0100 Subject: [PATCH 16/22] Update test outputs --- .../airy_isostasy_initial_topo/screen-output | 23 +- .../screen-output | 14 +- .../statistics | 2 +- .../screen-output | 17 +- .../solution/solution-00000.0000.gnuplot | 390 +----------------- .../statistics | 2 +- 6 files changed, 51 insertions(+), 397 deletions(-) diff --git a/tests/airy_isostasy_initial_topo/screen-output b/tests/airy_isostasy_initial_topo/screen-output index e708adb2a6a..e13824c064d 100644 --- a/tests/airy_isostasy_initial_topo/screen-output +++ b/tests/airy_isostasy_initial_topo/screen-output @@ -1,4 +1,6 @@ +Loading shared library <./libairy_isostasy_initial_topo.debug.so> + Number of active cells: 1,024 (on 6 levels) Number of degrees of freedom: 17,989 (8,450+1,089+4,225+4,225) @@ -38,9 +40,10 @@ Number of mesh deformation degrees of freedom: 4,162 Solving Stokes system... 200+16 iterations. Postprocessing: - RMS, max velocity: 1.76e+06 m/year, 2.3e+06 m/year - Pressure min/avg/max: -1.132 Pa, 5.012 Pa, 10 Pa - Topography min/max: 0 m, 0.05 m + RMS, max velocity: 1.76e+06 m/year, 2.3e+06 m/year + Pressure min/avg/max: -1.132 Pa, 5.012 Pa, 10 Pa + Topography min/max: 0 m, 0.05 m + Points lie in domain: truefalsefalse *** Timestep 1: t=1.69957e-09 years, dt=1.69957e-09 years Solving mesh displacement system... 8 iterations. @@ -50,9 +53,10 @@ Number of mesh deformation degrees of freedom: 4,162 Solving Stokes system... 200+10 iterations. Postprocessing: - RMS, max velocity: 1.62e+06 m/year, 2.11e+06 m/year - Pressure min/avg/max: -1.119 Pa, 5.004 Pa, 10.01 Pa - Topography min/max: -4.848e-05 m, 0.04609 m + RMS, max velocity: 1.62e+06 m/year, 2.11e+06 m/year + Pressure min/avg/max: -1.119 Pa, 5.004 Pa, 10.01 Pa + Topography min/max: -4.848e-05 m, 0.04609 m + Points lie in domain: truefalsefalse *** Timestep 2: t=3.5e-09 years, dt=1.80043e-09 years Solving mesh displacement system... 8 iterations. @@ -62,9 +66,10 @@ Number of mesh deformation degrees of freedom: 4,162 Solving Stokes system... 200+4 iterations. Postprocessing: - RMS, max velocity: 1.53e+06 m/year, 1.99e+06 m/year - Pressure min/avg/max: -0.8445 Pa, 5.003 Pa, 10.02 Pa - Topography min/max: -9.971e-05 m, 0.04229 m + RMS, max velocity: 1.53e+06 m/year, 1.99e+06 m/year + Pressure min/avg/max: -0.8445 Pa, 5.003 Pa, 10.02 Pa + Topography min/max: -9.971e-05 m, 0.04229 m + Points lie in domain: truefalsefalse Termination requested by criterion: end time diff --git a/tests/airy_isostasy_initial_topo_lithostatic_pressure/screen-output b/tests/airy_isostasy_initial_topo_lithostatic_pressure/screen-output index e5d91664954..4d894d2c8ba 100644 --- a/tests/airy_isostasy_initial_topo_lithostatic_pressure/screen-output +++ b/tests/airy_isostasy_initial_topo_lithostatic_pressure/screen-output @@ -9,7 +9,7 @@ Number of mesh deformation degrees of freedom: 2,178 Skipping temperature solve because RHS is zero. Solving C_1 system ... 0 iterations. Rebuilding Stokes preconditioner... - Solving Stokes system... 200+8 iterations. + Solving Stokes system... 200+5 iterations. Number of active cells: 1,204 (on 7 levels) Number of degrees of freedom: 21,360 (10,034+1,292+5,017+5,017) @@ -21,24 +21,26 @@ Number of mesh deformation degrees of freedom: 2,584 Skipping temperature solve because RHS is zero. Solving C_1 system ... 0 iterations. Rebuilding Stokes preconditioner... - Solving Stokes system... 167+0 iterations. + Solving Stokes system... 169+0 iterations. Number of active cells: 1,948 (on 8 levels) Number of degrees of freedom: 34,621 (16,270+2,081+8,135+8,135) Number of mesh deformation degrees of freedom: 4,162 Solving mesh displacement system... 0 iterations. + Added initial topography to grid + *** Timestep 0: t=0 years, dt=0 years Solving mesh displacement system... 0 iterations. Skipping temperature solve because RHS is zero. Solving C_1 system ... 0 iterations. Rebuilding Stokes preconditioner... - Solving Stokes system... 200+6 iterations. + Solving Stokes system... 200+16 iterations. Postprocessing: - RMS, max velocity: 6.63e+05 m/year, 8.88e+05 m/year - Pressure min/avg/max: -0.6507 Pa, 5.005 Pa, 10 Pa - Topography min/max: 0 m, 0 m + RMS, max velocity: 9.35e+05 m/year, 1.25e+06 m/year + Pressure min/avg/max: -0.5111 Pa, 5.145 Pa, 10.29 Pa + Topography min/max: 0 m, 0.05 m Termination requested by criterion: end time diff --git a/tests/airy_isostasy_initial_topo_lithostatic_pressure/statistics b/tests/airy_isostasy_initial_topo_lithostatic_pressure/statistics index 70cc69bf902..c9136bd4bbc 100644 --- a/tests/airy_isostasy_initial_topo_lithostatic_pressure/statistics +++ b/tests/airy_isostasy_initial_topo_lithostatic_pressure/statistics @@ -17,4 +17,4 @@ # 17: Maximal pressure (Pa) # 18: Minimum topography (m) # 19: Maximum topography (m) -0 0.000000000000e+00 0.000000000000e+00 1948 18351 8135 8135 0 0 206 444 840 6.62545995e+05 8.87572504e+05 -6.50724076e-01 5.00501899e+00 1.00000000e+01 0.00000000e+00 0.00000000e+00 +0 0.000000000000e+00 0.000000000000e+00 1948 18351 8135 8135 0 0 216 809 882 9.35492994e+05 1.25102336e+06 -5.11120571e-01 5.14534156e+00 1.02938926e+01 0.00000000e+00 5.00000000e-02 diff --git a/tests/chunk_lithostatic_pressure_2d_initial_topography/screen-output b/tests/chunk_lithostatic_pressure_2d_initial_topography/screen-output index 3fc22aca562..932104ac7c6 100644 --- a/tests/chunk_lithostatic_pressure_2d_initial_topography/screen-output +++ b/tests/chunk_lithostatic_pressure_2d_initial_topography/screen-output @@ -1,17 +1,20 @@ -Number of active cells: 8 (on 2 levels) -Number of degrees of freedom: 652 (375+27+125+125) + + Loading Ascii data boundary file ASPECT_DIR/data/geometry-model/initial-topography-model/ascii-data/test/shell_2d_outer.0.txt. + +Number of active cells: 4 (on 2 levels) +Number of degrees of freedom: 109 (50+9+25+25) *** Timestep 0: t=0 years, dt=0 years Solving temperature system... 0 iterations. - Solving C_1 system ... 11 iterations. + Solving C_1 system ... 8 iterations. Rebuilding Stokes preconditioner... - Solving Stokes system... 22+0 iterations. + Solving Stokes system... 19+0 iterations. Postprocessing: - RMS, max velocity: 0.309 m/year, 1.86 m/year - Pressure min/avg/max: -1.965e+07 Pa, 3.29e+10 Pa, 9.157e+10 Pa - Writing graphical output: output-chunk_lithostatic_pressure_2d/solution/solution-00000 + RMS, max velocity: 0.275 m/year, 0.428 m/year + Pressure min/avg/max: -2.14e+09 Pa, 3.928e+10 Pa, 9.249e+10 Pa + Writing graphical output: output-chunk_lithostatic_pressure_2d_initial_topography/solution/solution-00000 Termination requested by criterion: end time diff --git a/tests/chunk_lithostatic_pressure_2d_initial_topography/solution/solution-00000.0000.gnuplot b/tests/chunk_lithostatic_pressure_2d_initial_topography/solution/solution-00000.0000.gnuplot index baca9c2fcbc..66d4ba74b7d 100644 --- a/tests/chunk_lithostatic_pressure_2d_initial_topography/solution/solution-00000.0000.gnuplot +++ b/tests/chunk_lithostatic_pressure_2d_initial_topography/solution/solution-00000.0000.gnuplot @@ -4,388 +4,32 @@ # # For a description of the GNUPLOT format see the GNUPLOT manual. # -#

-3.41827e+06 602733 0 0 0 0 9.14242e+10 1600 0 9.14598e+10 4.13146e+06 3.14744e+07 4.13146e+06 9.14342e+10 5.54979e+06 3.14744e+07 5.54979e+06 9.10905e+10 1.04362e-13 4.50663e+07 -4.84624e+06 854523 0 -0.057618 -0.0101596 0 4.56474e+10 1600 -1.84884e-16 4.56173e+10 -6.04158e+06 8.6398e+06 -6.04158e+06 4.56496e+10 1.52343e+06 8.6398e+06 1.52343e+06 4.56317e+10 9.68689e-15 -1.78755e+07 +#

+3.471e+06 0 0 0 9.2485e+10 1600 0 9.24304e+10 -2.48143e+06 -2.48143e+06 9.24992e+10 1.72466e-14 1.10592e+09 +4.966e+06 0 0.66533 0 4.51592e+10 1600 -3.0288e-17 4.51574e+10 4.60252e+07 4.60252e+07 4.52082e+10 2.62907e-14 9.12587e+08 +3.14579e+06 1.46691e+06 -0.0828562 0.177686 9.24905e+10 1600 0 9.24651e+10 -351844 -351844 9.24576e+10 1.90575e-15 1.20296e+09 +4.47532e+06 2.08688e+06 -0.180914 0.503458 4.61686e+10 1600 -3.6284e-17 4.62145e+10 -1.54555e+06 -1.54555e+06 4.60916e+10 3.07426e-14 1.03836e+09 -3.41827e+06 602733 0 0 0 0 9.14242e+10 1600 0 9.14598e+10 4.13146e+06 3.14744e+07 4.13146e+06 9.14342e+10 5.54979e+06 3.14744e+07 5.54979e+06 9.10905e+10 1.04362e-13 4.50663e+07 -3.00597e+06 1.7355e+06 0 0.0355464 -0.0615682 0 9.14379e+10 1600 0 9.1461e+10 1.6499e+07 2.71656e+07 1.6499e+07 9.1447e+10 1.7322e+07 2.71656e+07 1.7322e+07 9.1101e+10 1.03565e-13 5.8845e+07 +4.966e+06 0 0.66533 0 4.51592e+10 1600 -3.0288e-17 4.51422e+10 4.59544e+07 4.59544e+07 4.52082e+10 2.82955e-14 9.12587e+08 +6.461e+06 0 0 0 -2.13979e+09 1600 1 7.12564e+10 4.62219e+09 4.62219e+09 -3.37012e+10 2.6341e-14 -2.13979e+09 -3.41827e+06 602733 0 0 0 0 9.14242e+10 1600 0 9.14598e+10 4.13146e+06 3.14744e+07 4.13146e+06 9.14342e+10 5.54979e+06 3.14744e+07 5.54979e+06 9.10905e+10 1.04362e-13 4.50663e+07 -3.21212e+06 566384 1.18715e+06 0.185945 0.0327872 -0.518762 9.13123e+10 1600 0 9.14457e+10 -3.67688e+07 -1.49744e+08 -3.67688e+07 9.12698e+10 1.57118e+08 -1.49744e+08 1.57118e+08 9.16547e+10 1.46277e-13 -6.68253e+07 +4.47532e+06 2.08688e+06 -0.180914 0.503458 4.61686e+10 1600 -3.6284e-17 4.61908e+10 -2.10842e+07 -2.10842e+07 4.60783e+10 3.00476e-14 1.03836e+09 +5.80485e+06 2.70685e+06 0.0146399 -0.0293751 -4.47459e+07 1600 1 -1.56193e+10 2.42465e+10 2.42465e+10 6.92338e+10 2.44331e-14 -4.47459e+07 -4.84624e+06 854523 0 -0.057618 -0.0101596 0 4.56474e+10 1600 -1.84884e-16 4.56173e+10 -6.04158e+06 8.6398e+06 -6.04158e+06 4.56496e+10 1.52343e+06 8.6398e+06 1.52343e+06 4.56317e+10 9.68689e-15 -1.78755e+07 -4.26171e+06 2.4605e+06 0 -0.000689281 -0.112987 0 4.56467e+10 1600 -9.13492e-17 4.56287e+10 -1.6118e+07 7.29339e+06 -1.6118e+07 4.56415e+10 4.68681e+06 7.29339e+06 4.68681e+06 4.56313e+10 9.75342e-15 -1.85507e+07 +3.14579e+06 1.46691e+06 -0.0828562 0.177686 9.24905e+10 1600 0 9.24814e+10 -3.06337e+07 -3.06337e+07 9.25126e+10 1.7192e-14 1.20296e+09 +4.47532e+06 2.08688e+06 -0.180914 0.503458 4.61686e+10 1600 -3.6284e-17 4.6216e+10 -2.98026e+07 -2.98026e+07 4.6206e+10 1.51079e-14 1.03836e+09 +2.23112e+06 2.65894e+06 -0.60096 0.504265 9.24572e+10 1600 0 9.23332e+10 3.56637e+07 3.56637e+07 9.24968e+10 4.46275e-14 1.07814e+09 +3.15254e+06 3.75706e+06 -0.208402 -0.424846 4.7355e+10 1600 0 4.73654e+10 3.56775e+07 3.56775e+07 4.73449e+10 1.85572e-14 1.1692e+09 -4.84624e+06 854523 0 -0.057618 -0.0101596 0 4.56474e+10 1600 -1.84884e-16 4.56173e+10 -6.04158e+06 8.6398e+06 -6.04158e+06 4.56496e+10 1.52343e+06 8.6398e+06 1.52343e+06 4.56317e+10 9.68689e-15 -1.78755e+07 -4.55397e+06 802989 1.68308e+06 -0.148579 -0.0261985 -0.247956 4.56703e+10 1600 -5.52646e-17 4.56203e+10 -1.9544e+07 -2.35491e+06 -1.9544e+07 4.56633e+10 -223785 -2.35491e+06 -223785 4.57163e+10 2.59646e-14 5.05121e+06 +4.47532e+06 2.08688e+06 -0.180914 0.503458 4.61686e+10 1600 -3.6284e-17 4.61923e+10 -4.93617e+07 -4.93617e+07 4.61926e+10 2.4681e-14 1.03836e+09 +5.80485e+06 2.70685e+06 0.0146399 -0.0293751 -4.47459e+07 1600 1 -2.82182e+10 4.97887e+10 4.97887e+10 1.76241e+10 2.74057e-14 -4.47459e+07 -3.00597e+06 1.7355e+06 0 0.0355464 -0.0615682 0 9.14379e+10 1600 0 9.1461e+10 1.6499e+07 2.71656e+07 1.6499e+07 9.1447e+10 1.7322e+07 2.71656e+07 1.7322e+07 9.1101e+10 1.03565e-13 5.8845e+07 -4.26171e+06 2.4605e+06 0 -0.000689281 -0.112987 0 4.56467e+10 1600 -9.13492e-17 4.56287e+10 -1.6118e+07 7.29339e+06 -1.6118e+07 4.56415e+10 4.68681e+06 7.29339e+06 4.68681e+06 4.56313e+10 9.75342e-15 -1.85507e+07 - - -3.00597e+06 1.7355e+06 0 0.0355464 -0.0615682 0 9.14379e+10 1600 0 9.1461e+10 1.6499e+07 2.71656e+07 1.6499e+07 9.1447e+10 1.7322e+07 2.71656e+07 1.7322e+07 9.1101e+10 1.03565e-13 5.8845e+07 -2.82469e+06 1.63084e+06 1.18715e+06 0.19034 0.0376886 -0.504666 9.13093e+10 1600 0 9.13328e+10 9.30997e+07 -1.30653e+07 9.30997e+07 9.13805e+10 -2.1834e+08 -1.30653e+07 -2.1834e+08 9.16534e+10 1.46993e-13 -6.9817e+07 - - -4.26171e+06 2.4605e+06 0 -0.000689281 -0.112987 0 4.56467e+10 1600 -9.13492e-17 4.56287e+10 -1.6118e+07 7.29339e+06 -1.6118e+07 4.56415e+10 4.68681e+06 7.29339e+06 4.68681e+06 4.56313e+10 9.75342e-15 -1.85507e+07 -4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56204e+10 -2.00601e+07 684454 -2.00601e+07 4.56656e+10 -4.28112e+06 684454 -4.28112e+06 4.57146e+10 2.56922e-14 4.8163e+06 - - -3.21212e+06 566384 1.18715e+06 0.185945 0.0327872 -0.518762 9.13123e+10 1600 0 9.14457e+10 -3.67688e+07 -1.49744e+08 -3.67688e+07 9.12698e+10 1.57118e+08 -1.49744e+08 1.57118e+08 9.16547e+10 1.46277e-13 -6.68253e+07 -4.55397e+06 802989 1.68308e+06 -0.148579 -0.0261985 -0.247956 4.56703e+10 1600 -5.52646e-17 4.56203e+10 -1.9544e+07 -2.35491e+06 -1.9544e+07 4.56633e+10 -223785 -2.35491e+06 -223785 4.57163e+10 2.59646e-14 5.05121e+06 - - -3.21212e+06 566384 1.18715e+06 0.185945 0.0327872 -0.518762 9.13123e+10 1600 0 9.14457e+10 -3.67688e+07 -1.49744e+08 -3.67688e+07 9.12698e+10 1.57118e+08 -1.49744e+08 1.57118e+08 9.16547e+10 1.46277e-13 -6.68253e+07 -2.82469e+06 1.63084e+06 1.18715e+06 0.19034 0.0376886 -0.504666 9.13093e+10 1600 0 9.13328e+10 9.30997e+07 -1.30653e+07 9.30997e+07 9.13805e+10 -2.1834e+08 -1.30653e+07 -2.1834e+08 9.16534e+10 1.46993e-13 -6.9817e+07 - - -4.55397e+06 802989 1.68308e+06 -0.148579 -0.0261985 -0.247956 4.56703e+10 1600 -5.52646e-17 4.56203e+10 -1.9544e+07 -2.35491e+06 -1.9544e+07 4.56633e+10 -223785 -2.35491e+06 -223785 4.57163e+10 2.59646e-14 5.05121e+06 -4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56204e+10 -2.00601e+07 684454 -2.00601e+07 4.56656e+10 -4.28112e+06 684454 -4.28112e+06 4.57146e+10 2.56922e-14 4.8163e+06 - - -2.82469e+06 1.63084e+06 1.18715e+06 0.19034 0.0376886 -0.504666 9.13093e+10 1600 0 9.13328e+10 9.30997e+07 -1.30653e+07 9.30997e+07 9.13805e+10 -2.1834e+08 -1.30653e+07 -2.1834e+08 9.16534e+10 1.46993e-13 -6.9817e+07 -4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56204e+10 -2.00601e+07 684454 -2.00601e+07 4.56656e+10 -4.28112e+06 684454 -4.28112e+06 4.57146e+10 2.56922e-14 4.8163e+06 - - -4.84624e+06 854523 0 -0.057618 -0.0101596 0 4.56474e+10 1600 -1.84884e-16 4.56422e+10 -1.65431e+06 8.63978e+06 -1.65431e+06 4.56504e+10 1.52343e+06 8.63978e+06 1.52343e+06 4.56317e+10 6.45996e-15 -1.78755e+07 -6.27421e+06 1.10631e+06 0 0 0 0 6.38156e+06 1600 1 4.80547e+08 -1.17027e+08 9.8169e+07 -1.17027e+08 6.31248e+08 1.73098e+07 9.8169e+07 1.73098e+07 3.1943e+08 1.09486e-16 6.38156e+06 - - -4.84624e+06 854523 0 -0.057618 -0.0101596 0 4.56474e+10 1600 -1.84884e-16 4.56422e+10 -1.65431e+06 8.63978e+06 -1.65431e+06 4.56504e+10 1.52343e+06 8.63978e+06 1.52343e+06 4.56317e+10 6.45996e-15 -1.78755e+07 -4.26171e+06 2.4605e+06 0 -0.000689281 -0.112987 0 4.56467e+10 1600 -9.13492e-17 4.56486e+10 -6.88188e+06 7.29337e+06 -6.88188e+06 4.56455e+10 4.68681e+06 7.29337e+06 4.68681e+06 4.56313e+10 7.1962e-15 -1.85507e+07 - - -4.84624e+06 854523 0 -0.057618 -0.0101596 0 4.56474e+10 1600 -1.84884e-16 4.56422e+10 -1.65431e+06 8.63978e+06 -1.65431e+06 4.56504e+10 1.52343e+06 8.63978e+06 1.52343e+06 4.56317e+10 6.45996e-15 -1.78755e+07 -4.55397e+06 802989 1.68308e+06 -0.148579 -0.0261985 -0.247956 4.56703e+10 1600 -5.52646e-17 4.56596e+10 -1.26263e+07 -2.35403e+07 -1.26263e+07 4.56645e+10 -3.95929e+06 -2.35403e+07 -3.95929e+06 4.56952e+10 1.66034e-14 5.05121e+06 - - -6.27421e+06 1.10631e+06 0 0 0 0 6.38156e+06 1600 1 4.80547e+08 -1.17027e+08 9.8169e+07 -1.17027e+08 6.31248e+08 1.73098e+07 9.8169e+07 1.73098e+07 3.1943e+08 1.09486e-16 6.38156e+06 -5.51745e+06 3.1855e+06 0 -0.00516999 0.00895468 0 1.11727e+07 1600 1 -1.34887e+08 3.51473e+08 6.55746e+07 3.51473e+08 -1.22425e+09 5.55838e+07 6.55746e+07 5.55838e+07 2.56716e+08 4.24209e-16 1.11727e+07 - - -6.27421e+06 1.10631e+06 0 0 0 0 6.38156e+06 1600 1 4.80547e+08 -1.17027e+08 9.8169e+07 -1.17027e+08 6.31248e+08 1.73098e+07 9.8169e+07 1.73098e+07 3.1943e+08 1.09486e-16 6.38156e+06 -5.89583e+06 1.03959e+06 2.17901e+06 0.0170632 0.00300871 -0.0476041 -3.4701e+06 1600 1 5.00929e+09 1.10033e+09 1.30345e+09 1.10033e+09 -1.16799e+09 1.40563e+08 1.30345e+09 1.40563e+08 3.27351e+09 1.80826e-15 -3.4701e+06 - - -4.26171e+06 2.4605e+06 0 -0.000689281 -0.112987 0 4.56467e+10 1600 -9.13492e-17 4.56486e+10 -6.88188e+06 7.29337e+06 -6.88188e+06 4.56455e+10 4.68681e+06 7.29337e+06 4.68681e+06 4.56313e+10 7.1962e-15 -1.85507e+07 -5.51745e+06 3.1855e+06 0 -0.00516999 0.00895468 0 1.11727e+07 1600 1 -1.34887e+08 3.51473e+08 6.55746e+07 3.51473e+08 -1.22425e+09 5.55838e+07 6.55746e+07 5.55838e+07 2.56716e+08 4.24209e-16 1.11727e+07 - - -4.26171e+06 2.4605e+06 0 -0.000689281 -0.112987 0 4.56467e+10 1600 -9.13492e-17 4.56486e+10 -6.88188e+06 7.29337e+06 -6.88188e+06 4.56455e+10 4.68681e+06 7.29337e+06 4.68681e+06 4.56313e+10 7.1962e-15 -1.85507e+07 -4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56509e+10 -3.6295e+06 -1.76123e+07 -3.6295e+06 4.56744e+10 -1.53349e+07 -1.76123e+07 -1.53349e+07 4.56938e+10 1.59777e-14 4.8163e+06 - - -5.51745e+06 3.1855e+06 0 -0.00516999 0.00895468 0 1.11727e+07 1600 1 -1.34887e+08 3.51473e+08 6.55746e+07 3.51473e+08 -1.22425e+09 5.55838e+07 6.55746e+07 5.55838e+07 2.56716e+08 4.24209e-16 1.11727e+07 -5.18471e+06 2.99339e+06 2.17901e+06 -0.0171684 0.0586576 -0.0397298 -1.83017e+06 1600 1 3.17549e+09 3.75887e+09 2.18939e+09 3.75887e+09 -3.25586e+08 -6.23188e+08 2.18939e+09 -6.23188e+08 2.94865e+09 2.40568e-15 -1.83017e+06 - - -4.55397e+06 802989 1.68308e+06 -0.148579 -0.0261985 -0.247956 4.56703e+10 1600 -5.52646e-17 4.56596e+10 -1.26263e+07 -2.35403e+07 -1.26263e+07 4.56645e+10 -3.95929e+06 -2.35403e+07 -3.95929e+06 4.56952e+10 1.66034e-14 5.05121e+06 -5.89583e+06 1.03959e+06 2.17901e+06 0.0170632 0.00300871 -0.0476041 -3.4701e+06 1600 1 5.00929e+09 1.10033e+09 1.30345e+09 1.10033e+09 -1.16799e+09 1.40563e+08 1.30345e+09 1.40563e+08 3.27351e+09 1.80826e-15 -3.4701e+06 - - -4.55397e+06 802989 1.68308e+06 -0.148579 -0.0261985 -0.247956 4.56703e+10 1600 -5.52646e-17 4.56596e+10 -1.26263e+07 -2.35403e+07 -1.26263e+07 4.56645e+10 -3.95929e+06 -2.35403e+07 -3.95929e+06 4.56952e+10 1.66034e-14 5.05121e+06 -4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56509e+10 -3.6295e+06 -1.76123e+07 -3.6295e+06 4.56744e+10 -1.53349e+07 -1.76123e+07 -1.53349e+07 4.56938e+10 1.59777e-14 4.8163e+06 - - -5.89583e+06 1.03959e+06 2.17901e+06 0.0170632 0.00300871 -0.0476041 -3.4701e+06 1600 1 5.00929e+09 1.10033e+09 1.30345e+09 1.10033e+09 -1.16799e+09 1.40563e+08 1.30345e+09 1.40563e+08 3.27351e+09 1.80826e-15 -3.4701e+06 -5.18471e+06 2.99339e+06 2.17901e+06 -0.0171684 0.0586576 -0.0397298 -1.83017e+06 1600 1 3.17549e+09 3.75887e+09 2.18939e+09 3.75887e+09 -3.25586e+08 -6.23188e+08 2.18939e+09 -6.23188e+08 2.94865e+09 2.40568e-15 -1.83017e+06 - - -4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56509e+10 -3.6295e+06 -1.76123e+07 -3.6295e+06 4.56744e+10 -1.53349e+07 -1.76123e+07 -1.53349e+07 4.56938e+10 1.59777e-14 4.8163e+06 -5.18471e+06 2.99339e+06 2.17901e+06 -0.0171684 0.0586576 -0.0397298 -1.83017e+06 1600 1 3.17549e+09 3.75887e+09 2.18939e+09 3.75887e+09 -3.25586e+08 -6.23188e+08 2.18939e+09 -6.23188e+08 2.94865e+09 2.40568e-15 -1.83017e+06 - - -3.00597e+06 1.7355e+06 0 0.0355464 -0.0615682 0 9.14379e+10 1600 0 9.1468e+10 7.93689e+06 2.71656e+07 7.93689e+06 9.14555e+10 1.7322e+07 2.71656e+07 1.7322e+07 9.1101e+10 1.05519e-13 5.8845e+07 -4.26171e+06 2.4605e+06 0 -0.000689281 -0.112987 0 4.56467e+10 1600 -9.13492e-17 4.56324e+10 -1.8969e+07 7.29339e+06 -1.8969e+07 4.564e+10 4.68681e+06 7.29339e+06 4.68681e+06 4.56313e+10 1.06925e-14 -1.85507e+07 - - -3.00597e+06 1.7355e+06 0 0.0355464 -0.0615682 0 9.14379e+10 1600 0 9.1468e+10 7.93689e+06 2.71656e+07 7.93689e+06 9.14555e+10 1.7322e+07 2.71656e+07 1.7322e+07 9.1101e+10 1.05519e-13 5.8845e+07 -2.23112e+06 2.65894e+06 0 -0.425883 0.357358 0 9.13978e+10 1600 0 9.13185e+10 5.12101e+07 761904 5.12101e+07 9.14355e+10 4.14712e+07 761904 4.14712e+07 9.10642e+10 1.00487e-13 1.87398e+07 - - -3.00597e+06 1.7355e+06 0 0.0355464 -0.0615682 0 9.14379e+10 1600 0 9.1468e+10 7.93689e+06 2.71656e+07 7.93689e+06 9.14555e+10 1.7322e+07 2.71656e+07 1.7322e+07 9.1101e+10 1.05519e-13 5.8845e+07 -2.82469e+06 1.63084e+06 1.18715e+06 0.19034 0.0376886 -0.504666 9.13093e+10 1600 0 9.14467e+10 2.54924e+07 -1.92175e+08 2.54924e+07 9.12729e+10 9.18863e+07 -1.92175e+08 9.18863e+07 9.16534e+10 1.43442e-13 -6.9817e+07 - - -4.26171e+06 2.4605e+06 0 -0.000689281 -0.112987 0 4.56467e+10 1600 -9.13492e-17 4.56324e+10 -1.8969e+07 7.29339e+06 -1.8969e+07 4.564e+10 4.68681e+06 7.29339e+06 4.68681e+06 4.56313e+10 1.06925e-14 -1.85507e+07 -3.16316e+06 3.7697e+06 0 0.0866503 -0.288584 0 4.56761e+10 1600 0 4.56746e+10 -1.27688e+07 5.09732e+06 -1.27688e+07 4.56554e+10 6.93062e+06 5.09732e+06 6.93062e+06 4.56659e+10 9.07151e-15 1.08148e+07 - - -4.26171e+06 2.4605e+06 0 -0.000689281 -0.112987 0 4.56467e+10 1600 -9.13492e-17 4.56324e+10 -1.8969e+07 7.29339e+06 -1.8969e+07 4.564e+10 4.68681e+06 7.29339e+06 4.68681e+06 4.56313e+10 1.06925e-14 -1.85507e+07 -4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56466e+10 -3.70481e+07 989103 -3.70481e+07 4.56457e+10 -4.80879e+06 989103 -4.80879e+06 4.57146e+10 2.71892e-14 4.8163e+06 - - -2.23112e+06 2.65894e+06 0 -0.425883 0.357358 0 9.13978e+10 1600 0 9.13185e+10 5.12101e+07 761904 5.12101e+07 9.14355e+10 4.14712e+07 761904 4.14712e+07 9.10642e+10 1.00487e-13 1.87398e+07 -3.16316e+06 3.7697e+06 0 0.0866503 -0.288584 0 4.56761e+10 1600 0 4.56746e+10 -1.27688e+07 5.09732e+06 -1.27688e+07 4.56554e+10 6.93062e+06 5.09732e+06 6.93062e+06 4.56659e+10 9.07151e-15 1.08148e+07 - - -2.23112e+06 2.65894e+06 0 -0.425883 0.357358 0 9.13978e+10 1600 0 9.13185e+10 5.12101e+07 761904 5.12101e+07 9.14355e+10 4.14712e+07 761904 4.14712e+07 9.10642e+10 1.00487e-13 1.87398e+07 -2.09656e+06 2.49859e+06 1.18715e+06 0.151969 0.147411 -0.578639 9.13366e+10 1600 0 9.12783e+10 6.4997e+07 4.33178e+07 6.4997e+07 9.14735e+10 -1.96226e+08 4.33178e+07 -1.96226e+08 9.16845e+10 1.46516e-13 -4.25455e+07 - - -3.16316e+06 3.7697e+06 0 0.0866503 -0.288584 0 4.56761e+10 1600 0 4.56746e+10 -1.27688e+07 5.09732e+06 -1.27688e+07 4.56554e+10 6.93062e+06 5.09732e+06 6.93062e+06 4.56659e+10 9.07151e-15 1.08148e+07 -2.9724e+06 3.54236e+06 1.68308e+06 -0.120017 -0.203298 -0.283249 4.56924e+10 1600 0 4.56678e+10 -1.6427e+07 8.96147e+06 -1.6427e+07 4.56519e+10 -1.03223e+07 8.96147e+06 -1.03223e+07 4.57353e+10 2.45702e-14 2.71349e+07 - - -2.82469e+06 1.63084e+06 1.18715e+06 0.19034 0.0376886 -0.504666 9.13093e+10 1600 0 9.14467e+10 2.54924e+07 -1.92175e+08 2.54924e+07 9.12729e+10 9.18863e+07 -1.92175e+08 9.18863e+07 9.16534e+10 1.43442e-13 -6.9817e+07 -4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56466e+10 -3.70481e+07 989103 -3.70481e+07 4.56457e+10 -4.80879e+06 989103 -4.80879e+06 4.57146e+10 2.71892e-14 4.8163e+06 - - -2.82469e+06 1.63084e+06 1.18715e+06 0.19034 0.0376886 -0.504666 9.13093e+10 1600 0 9.14467e+10 2.54924e+07 -1.92175e+08 2.54924e+07 9.12729e+10 9.18863e+07 -1.92175e+08 9.18863e+07 9.16534e+10 1.43442e-13 -6.9817e+07 -2.09656e+06 2.49859e+06 1.18715e+06 0.151969 0.147411 -0.578639 9.13366e+10 1600 0 9.12783e+10 6.4997e+07 4.33178e+07 6.4997e+07 9.14735e+10 -1.96226e+08 4.33178e+07 -1.96226e+08 9.16845e+10 1.46516e-13 -4.25455e+07 - - -4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56466e+10 -3.70481e+07 989103 -3.70481e+07 4.56457e+10 -4.80879e+06 989103 -4.80879e+06 4.57146e+10 2.71892e-14 4.8163e+06 -2.9724e+06 3.54236e+06 1.68308e+06 -0.120017 -0.203298 -0.283249 4.56924e+10 1600 0 4.56678e+10 -1.6427e+07 8.96147e+06 -1.6427e+07 4.56519e+10 -1.03223e+07 8.96147e+06 -1.03223e+07 4.57353e+10 2.45702e-14 2.71349e+07 - - -2.09656e+06 2.49859e+06 1.18715e+06 0.151969 0.147411 -0.578639 9.13366e+10 1600 0 9.12783e+10 6.4997e+07 4.33178e+07 6.4997e+07 9.14735e+10 -1.96226e+08 4.33178e+07 -1.96226e+08 9.16845e+10 1.46516e-13 -4.25455e+07 -2.9724e+06 3.54236e+06 1.68308e+06 -0.120017 -0.203298 -0.283249 4.56924e+10 1600 0 4.56678e+10 -1.6427e+07 8.96147e+06 -1.6427e+07 4.56519e+10 -1.03223e+07 8.96147e+06 -1.03223e+07 4.57353e+10 2.45702e-14 2.71349e+07 - - -4.26171e+06 2.4605e+06 0 -0.000689281 -0.112987 0 4.56467e+10 1600 -9.13492e-17 4.56524e+10 -9.73295e+06 7.29337e+06 -9.73295e+06 4.5644e+10 4.68681e+06 7.29337e+06 4.68681e+06 4.56313e+10 8.39786e-15 -1.85507e+07 -5.51745e+06 3.1855e+06 0 -0.00516999 0.00895468 0 1.11727e+07 1600 1 -3.23725e+07 1.88972e+08 6.55746e+07 1.88972e+08 -9.6888e+08 5.55838e+07 6.55746e+07 5.55838e+07 2.56716e+08 3.36731e-16 1.11727e+07 - - -4.26171e+06 2.4605e+06 0 -0.000689281 -0.112987 0 4.56467e+10 1600 -9.13492e-17 4.56524e+10 -9.73295e+06 7.29337e+06 -9.73295e+06 4.5644e+10 4.68681e+06 7.29337e+06 4.68681e+06 4.56313e+10 8.39786e-15 -1.85507e+07 -3.16316e+06 3.7697e+06 0 0.0866503 -0.288584 0 4.56761e+10 1600 0 4.56788e+10 -8.8847e+06 5.09731e+06 -8.8847e+06 4.56587e+10 6.93062e+06 5.09731e+06 6.93062e+06 4.56659e+10 8.0177e-15 1.08148e+07 - - -4.26171e+06 2.4605e+06 0 -0.000689281 -0.112987 0 4.56467e+10 1600 -9.13492e-17 4.56524e+10 -9.73295e+06 7.29337e+06 -9.73295e+06 4.5644e+10 4.68681e+06 7.29337e+06 4.68681e+06 4.56313e+10 8.39786e-15 -1.85507e+07 -4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56771e+10 -2.06175e+07 -1.73077e+07 -2.06175e+07 4.56545e+10 -1.58625e+07 -1.73077e+07 -1.58625e+07 4.56938e+10 1.84777e-14 4.8163e+06 - - -5.51745e+06 3.1855e+06 0 -0.00516999 0.00895468 0 1.11727e+07 1600 1 -3.23725e+07 1.88972e+08 6.55746e+07 1.88972e+08 -9.6888e+08 5.55838e+07 6.55746e+07 5.55838e+07 2.56716e+08 3.36731e-16 1.11727e+07 -4.0952e+06 4.88047e+06 0 -0.063111 0.0529564 0 -1.38844e+07 1600 1 1.23456e+07 3.39785e+09 -2.70267e+08 3.39785e+09 1.07996e+07 2.20761e+08 -2.70267e+08 2.20761e+08 2.83632e+09 1.89255e-15 -1.38844e+07 - - -5.51745e+06 3.1855e+06 0 -0.00516999 0.00895468 0 1.11727e+07 1600 1 -3.23725e+07 1.88972e+08 6.55746e+07 1.88972e+08 -9.6888e+08 5.55838e+07 6.55746e+07 5.55838e+07 2.56716e+08 3.36731e-16 1.11727e+07 -5.18471e+06 2.99339e+06 2.17901e+06 -0.0171684 0.0586576 -0.0397298 -1.83017e+06 1600 1 3.28783e+09 3.56038e+09 2.14447e+09 3.56038e+09 2.49725e+07 -5.45396e+08 2.14447e+09 -5.45396e+08 2.94865e+09 2.27984e-15 -1.83017e+06 - - -3.16316e+06 3.7697e+06 0 0.0866503 -0.288584 0 4.56761e+10 1600 0 4.56788e+10 -8.8847e+06 5.09731e+06 -8.8847e+06 4.56587e+10 6.93062e+06 5.09731e+06 6.93062e+06 4.56659e+10 8.0177e-15 1.08148e+07 -4.0952e+06 4.88047e+06 0 -0.063111 0.0529564 0 -1.38844e+07 1600 1 1.23456e+07 3.39785e+09 -2.70267e+08 3.39785e+09 1.07996e+07 2.20761e+08 -2.70267e+08 2.20761e+08 2.83632e+09 1.89255e-15 -1.38844e+07 - - -3.16316e+06 3.7697e+06 0 0.0866503 -0.288584 0 4.56761e+10 1600 0 4.56788e+10 -8.8847e+06 5.09731e+06 -8.8847e+06 4.56587e+10 6.93062e+06 5.09731e+06 6.93062e+06 4.56659e+10 8.0177e-15 1.08148e+07 -2.9724e+06 3.54236e+06 1.68308e+06 -0.120017 -0.203298 -0.283249 4.56924e+10 1600 0 4.56611e+10 -7.03291e+06 -1.16048e+07 -7.03291e+06 4.56839e+10 -2.49822e+07 -1.16048e+07 -2.49822e+07 4.57141e+10 1.94695e-14 2.71349e+07 - - -4.0952e+06 4.88047e+06 0 -0.063111 0.0529564 0 -1.38844e+07 1600 1 1.23456e+07 3.39785e+09 -2.70267e+08 3.39785e+09 1.07996e+07 2.20761e+08 -2.70267e+08 2.20761e+08 2.83632e+09 1.89255e-15 -1.38844e+07 -3.84823e+06 4.58614e+06 2.17901e+06 -0.0834525 0.113859 -0.0922569 -1.96489e+07 1600 1 3.3168e+08 5.08411e+09 1.54245e+09 5.08411e+09 4.54904e+09 1.87102e+09 1.54245e+09 1.87102e+09 4.00823e+09 3.04113e-15 -1.96489e+07 - - -4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56771e+10 -2.06175e+07 -1.73077e+07 -2.06175e+07 4.56545e+10 -1.58625e+07 -1.73077e+07 -1.58625e+07 4.56938e+10 1.84777e-14 4.8163e+06 -5.18471e+06 2.99339e+06 2.17901e+06 -0.0171684 0.0586576 -0.0397298 -1.83017e+06 1600 1 3.28783e+09 3.56038e+09 2.14447e+09 3.56038e+09 2.49725e+07 -5.45396e+08 2.14447e+09 -5.45396e+08 2.94865e+09 2.27984e-15 -1.83017e+06 - - -4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56771e+10 -2.06175e+07 -1.73077e+07 -2.06175e+07 4.56545e+10 -1.58625e+07 -1.73077e+07 -1.58625e+07 4.56938e+10 1.84777e-14 4.8163e+06 -2.9724e+06 3.54236e+06 1.68308e+06 -0.120017 -0.203298 -0.283249 4.56924e+10 1600 0 4.56611e+10 -7.03291e+06 -1.16048e+07 -7.03291e+06 4.56839e+10 -2.49822e+07 -1.16048e+07 -2.49822e+07 4.57141e+10 1.94695e-14 2.71349e+07 - - -5.18471e+06 2.99339e+06 2.17901e+06 -0.0171684 0.0586576 -0.0397298 -1.83017e+06 1600 1 3.28783e+09 3.56038e+09 2.14447e+09 3.56038e+09 2.49725e+07 -5.45396e+08 2.14447e+09 -5.45396e+08 2.94865e+09 2.27984e-15 -1.83017e+06 -3.84823e+06 4.58614e+06 2.17901e+06 -0.0834525 0.113859 -0.0922569 -1.96489e+07 1600 1 3.3168e+08 5.08411e+09 1.54245e+09 5.08411e+09 4.54904e+09 1.87102e+09 1.54245e+09 1.87102e+09 4.00823e+09 3.04113e-15 -1.96489e+07 - - -2.9724e+06 3.54236e+06 1.68308e+06 -0.120017 -0.203298 -0.283249 4.56924e+10 1600 0 4.56611e+10 -7.03291e+06 -1.16048e+07 -7.03291e+06 4.56839e+10 -2.49822e+07 -1.16048e+07 -2.49822e+07 4.57141e+10 1.94695e-14 2.71349e+07 -3.84823e+06 4.58614e+06 2.17901e+06 -0.0834525 0.113859 -0.0922569 -1.96489e+07 1600 1 3.3168e+08 5.08411e+09 1.54245e+09 5.08411e+09 4.54904e+09 1.87102e+09 1.54245e+09 1.87102e+09 4.00823e+09 3.04113e-15 -1.96489e+07 - - -3.21212e+06 566384 1.18715e+06 0.185945 0.0327872 -0.518762 9.13123e+10 1600 0 9.12813e+10 -6.57581e+07 2.70157e+08 -6.57581e+07 9.12647e+10 2.31157e+08 2.70157e+08 2.31157e+08 9.05914e+10 2.67246e-13 -6.68253e+07 -4.55397e+06 802989 1.68308e+06 -0.148579 -0.0261985 -0.247956 4.56703e+10 1600 -5.52646e-17 4.56345e+10 -1.70458e+07 -1.33254e+06 -1.70458e+07 4.56638e+10 -43540.7 -1.33254e+06 -43540.7 4.56003e+10 1.804e-14 5.05121e+06 - - -3.21212e+06 566384 1.18715e+06 0.185945 0.0327872 -0.518762 9.13123e+10 1600 0 9.12813e+10 -6.57581e+07 2.70157e+08 -6.57581e+07 9.12647e+10 2.31157e+08 2.70157e+08 2.31157e+08 9.05914e+10 2.67246e-13 -6.68253e+07 -2.82469e+06 1.63084e+06 1.18715e+06 0.19034 0.0376886 -0.504666 9.13093e+10 1600 0 9.12035e+10 2.02277e+07 3.60459e+08 2.02277e+07 9.13395e+10 -8.36812e+06 3.60459e+08 -8.36812e+06 9.05849e+10 2.70251e-13 -6.9817e+07 - - -3.21212e+06 566384 1.18715e+06 0.185945 0.0327872 -0.518762 9.13123e+10 1600 0 9.12813e+10 -6.57581e+07 2.70157e+08 -6.57581e+07 9.12647e+10 2.31157e+08 2.70157e+08 2.31157e+08 9.05914e+10 2.67246e-13 -6.68253e+07 -2.61854e+06 461720 2.23112e+06 0 0 0 9.15676e+10 1600 0 9.18492e+10 4.86143e+07 -3.34847e+08 4.86143e+07 9.15779e+10 -5.90422e+07 -3.34847e+08 -5.90422e+07 9.21224e+10 2.19142e-13 1.88529e+08 - - -4.55397e+06 802989 1.68308e+06 -0.148579 -0.0261985 -0.247956 4.56703e+10 1600 -5.52646e-17 4.56345e+10 -1.70458e+07 -1.33254e+06 -1.70458e+07 4.56638e+10 -43540.7 -1.33254e+06 -43540.7 4.56003e+10 1.804e-14 5.05121e+06 -4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56315e+10 -1.35334e+07 1.17341e+06 -1.35334e+07 4.56694e+10 -4.32639e+06 1.17341e+06 -4.32639e+06 4.55995e+10 1.89028e-14 4.8163e+06 - - -4.55397e+06 802989 1.68308e+06 -0.148579 -0.0261985 -0.247956 4.56703e+10 1600 -5.52646e-17 4.56345e+10 -1.70458e+07 -1.33254e+06 -1.70458e+07 4.56638e+10 -43540.7 -1.33254e+06 -43540.7 4.56003e+10 1.804e-14 5.05121e+06 -3.71243e+06 654602 3.16316e+06 -0.708369 -0.124905 -0.603561 4.56608e+10 1600 -1.94027e-17 4.56527e+10 -1.19489e+08 2.62342e+07 -1.19489e+08 4.55834e+10 -1.0272e+08 2.62342e+07 -1.0272e+08 4.57566e+10 9.09935e-14 -4.52843e+06 - - -2.82469e+06 1.63084e+06 1.18715e+06 0.19034 0.0376886 -0.504666 9.13093e+10 1600 0 9.12035e+10 2.02277e+07 3.60459e+08 2.02277e+07 9.13395e+10 -8.36812e+06 3.60459e+08 -8.36812e+06 9.05849e+10 2.70251e-13 -6.9817e+07 -4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56315e+10 -1.35334e+07 1.17341e+06 -1.35334e+07 4.56694e+10 -4.32639e+06 1.17341e+06 -4.32639e+06 4.55995e+10 1.89028e-14 4.8163e+06 - - -2.82469e+06 1.63084e+06 1.18715e+06 0.19034 0.0376886 -0.504666 9.13093e+10 1600 0 9.12035e+10 2.02277e+07 3.60459e+08 2.02277e+07 9.13395e+10 -8.36812e+06 3.60459e+08 -8.36812e+06 9.05849e+10 2.70251e-13 -6.9817e+07 -2.30271e+06 1.32947e+06 2.23112e+06 -0.045879 0.0794648 0 9.15743e+10 1600 0 9.17898e+10 1.34251e+08 -2.95352e+08 1.34251e+08 9.16379e+10 -1.71012e+08 -2.95352e+08 -1.71012e+08 9.21342e+10 2.23142e-13 1.95234e+08 - - -4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56315e+10 -1.35334e+07 1.17341e+06 -1.35334e+07 4.56694e+10 -4.32639e+06 1.17341e+06 -4.32639e+06 4.55995e+10 1.89028e-14 4.8163e+06 -3.26466e+06 1.88485e+06 3.16316e+06 -0.659004 -0.297874 -0.603859 4.56587e+10 1600 -1.34824e-18 4.55015e+10 5.59097e+07 -4.75938e+07 5.59097e+07 4.57279e+10 9.34442e+07 -4.75938e+07 9.34442e+07 4.57547e+10 9.14777e-14 -6.62725e+06 - - -2.61854e+06 461720 2.23112e+06 0 0 0 9.15676e+10 1600 0 9.18492e+10 4.86143e+07 -3.34847e+08 4.86143e+07 9.15779e+10 -5.90422e+07 -3.34847e+08 -5.90422e+07 9.21224e+10 2.19142e-13 1.88529e+08 -3.71243e+06 654602 3.16316e+06 -0.708369 -0.124905 -0.603561 4.56608e+10 1600 -1.94027e-17 4.56527e+10 -1.19489e+08 2.62342e+07 -1.19489e+08 4.55834e+10 -1.0272e+08 2.62342e+07 -1.0272e+08 4.57566e+10 9.09935e-14 -4.52843e+06 - - -2.61854e+06 461720 2.23112e+06 0 0 0 9.15676e+10 1600 0 9.18492e+10 4.86143e+07 -3.34847e+08 4.86143e+07 9.15779e+10 -5.90422e+07 -3.34847e+08 -5.90422e+07 9.21224e+10 2.19142e-13 1.88529e+08 -2.30271e+06 1.32947e+06 2.23112e+06 -0.045879 0.0794648 0 9.15743e+10 1600 0 9.17898e+10 1.34251e+08 -2.95352e+08 1.34251e+08 9.16379e+10 -1.71012e+08 -2.95352e+08 -1.71012e+08 9.21342e+10 2.23142e-13 1.95234e+08 - - -3.71243e+06 654602 3.16316e+06 -0.708369 -0.124905 -0.603561 4.56608e+10 1600 -1.94027e-17 4.56527e+10 -1.19489e+08 2.62342e+07 -1.19489e+08 4.55834e+10 -1.0272e+08 2.62342e+07 -1.0272e+08 4.57566e+10 9.09935e-14 -4.52843e+06 -3.26466e+06 1.88485e+06 3.16316e+06 -0.659004 -0.297874 -0.603859 4.56587e+10 1600 -1.34824e-18 4.55015e+10 5.59097e+07 -4.75938e+07 5.59097e+07 4.57279e+10 9.34442e+07 -4.75938e+07 9.34442e+07 4.57547e+10 9.14777e-14 -6.62725e+06 - - -2.30271e+06 1.32947e+06 2.23112e+06 -0.045879 0.0794648 0 9.15743e+10 1600 0 9.17898e+10 1.34251e+08 -2.95352e+08 1.34251e+08 9.16379e+10 -1.71012e+08 -2.95352e+08 -1.71012e+08 9.21342e+10 2.23142e-13 1.95234e+08 -3.26466e+06 1.88485e+06 3.16316e+06 -0.659004 -0.297874 -0.603859 4.56587e+10 1600 -1.34824e-18 4.55015e+10 5.59097e+07 -4.75938e+07 5.59097e+07 4.57279e+10 9.34442e+07 -4.75938e+07 9.34442e+07 4.57547e+10 9.14777e-14 -6.62725e+06 - - -4.55397e+06 802989 1.68308e+06 -0.148579 -0.0261985 -0.247956 4.56703e+10 1600 -5.52646e-17 4.56737e+10 -1.01281e+07 -2.2518e+07 -1.01281e+07 4.5665e+10 -3.77906e+06 -2.2518e+07 -3.77906e+06 4.55793e+10 2.89328e-14 5.05121e+06 -5.89583e+06 1.03959e+06 2.17901e+06 0.0170632 0.00300871 -0.0476041 -3.4701e+06 1600 1 4.58898e+09 1.02622e+09 2.57071e+09 1.02622e+09 -1.18106e+09 3.64013e+08 2.57071e+09 3.64013e+08 -5.26011e+08 2.10803e-15 -3.4701e+06 - - -4.55397e+06 802989 1.68308e+06 -0.148579 -0.0261985 -0.247956 4.56703e+10 1600 -5.52646e-17 4.56737e+10 -1.01281e+07 -2.2518e+07 -1.01281e+07 4.5665e+10 -3.77906e+06 -2.2518e+07 -3.77906e+06 4.55793e+10 2.89328e-14 5.05121e+06 -4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.5662e+10 2.89722e+06 -1.71234e+07 2.89722e+06 4.56783e+10 -1.53802e+07 -1.71234e+07 -1.53802e+07 4.55787e+10 2.91045e-14 4.8163e+06 - - -4.55397e+06 802989 1.68308e+06 -0.148579 -0.0261985 -0.247956 4.56703e+10 1600 -5.52646e-17 4.56737e+10 -1.01281e+07 -2.2518e+07 -1.01281e+07 4.5665e+10 -3.77906e+06 -2.2518e+07 -3.77906e+06 4.55793e+10 2.89328e-14 5.05121e+06 -3.71243e+06 654602 3.16316e+06 -0.708369 -0.124905 -0.603561 4.56608e+10 1600 -1.94027e-17 4.5615e+10 -1.26129e+08 -5.85141e+06 -1.26129e+08 4.55822e+10 -1.08378e+08 -5.85141e+06 -1.08378e+08 4.57293e+10 9.17168e-14 -4.52843e+06 - - -5.89583e+06 1.03959e+06 2.17901e+06 0.0170632 0.00300871 -0.0476041 -3.4701e+06 1600 1 4.58898e+09 1.02622e+09 2.57071e+09 1.02622e+09 -1.18106e+09 3.64013e+08 2.57071e+09 3.64013e+08 -5.26011e+08 2.10803e-15 -3.4701e+06 -5.18471e+06 2.99339e+06 2.17901e+06 -0.0171684 0.0586576 -0.0397298 -1.83017e+06 1600 1 2.74028e+09 3.63106e+09 3.43815e+09 3.63106e+09 -3.28105e+08 -2.93871e+08 3.43815e+09 -2.93871e+08 -5.94477e+08 2.6705e-15 -1.83017e+06 - - -5.89583e+06 1.03959e+06 2.17901e+06 0.0170632 0.00300871 -0.0476041 -3.4701e+06 1600 1 4.58898e+09 1.02622e+09 2.57071e+09 1.02622e+09 -1.18106e+09 3.64013e+08 2.57071e+09 3.64013e+08 -5.26011e+08 2.10803e-15 -3.4701e+06 -4.80632e+06 847485 4.0952e+06 0 0 0 1.81123e+07 1600 1 -1.96764e+10 -2.92287e+09 -1.55519e+10 -2.92287e+09 -4.19079e+09 -2.74221e+09 -1.55519e+10 -2.74221e+09 -1.46328e+10 8.94851e-15 1.81123e+07 - - -4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.5662e+10 2.89722e+06 -1.71234e+07 2.89722e+06 4.56783e+10 -1.53802e+07 -1.71234e+07 -1.53802e+07 4.55787e+10 2.91045e-14 4.8163e+06 -5.18471e+06 2.99339e+06 2.17901e+06 -0.0171684 0.0586576 -0.0397298 -1.83017e+06 1600 1 2.74028e+09 3.63106e+09 3.43815e+09 3.63106e+09 -3.28105e+08 -2.93871e+08 3.43815e+09 -2.93871e+08 -5.94477e+08 2.6705e-15 -1.83017e+06 - - -4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.5662e+10 2.89722e+06 -1.71234e+07 2.89722e+06 4.56783e+10 -1.53802e+07 -1.71234e+07 -1.53802e+07 4.55787e+10 2.91045e-14 4.8163e+06 -3.26466e+06 1.88485e+06 3.16316e+06 -0.659004 -0.297874 -0.603859 4.56587e+10 1600 -1.34824e-18 4.54752e+10 3.89659e+07 -7.38383e+07 3.89659e+07 4.57171e+10 7.66056e+07 -7.38383e+07 7.66056e+07 4.57286e+10 9.127e-14 -6.62725e+06 - - -5.18471e+06 2.99339e+06 2.17901e+06 -0.0171684 0.0586576 -0.0397298 -1.83017e+06 1600 1 2.74028e+09 3.63106e+09 3.43815e+09 3.63106e+09 -3.28105e+08 -2.93871e+08 3.43815e+09 -2.93871e+08 -5.94477e+08 2.6705e-15 -1.83017e+06 -4.22661e+06 2.44023e+06 4.0952e+06 -0.0644111 0.111563 3.93309e-17 1.73775e+07 1600 1 -1.68018e+10 -5.75744e+09 -1.36291e+10 -5.75744e+09 -9.28812e+09 -8.51327e+09 -1.36291e+10 -8.51327e+09 -1.41625e+10 8.74508e-15 1.73775e+07 - - -3.71243e+06 654602 3.16316e+06 -0.708369 -0.124905 -0.603561 4.56608e+10 1600 -1.94027e-17 4.5615e+10 -1.26129e+08 -5.85141e+06 -1.26129e+08 4.55822e+10 -1.08378e+08 -5.85141e+06 -1.08378e+08 4.57293e+10 9.17168e-14 -4.52843e+06 -4.80632e+06 847485 4.0952e+06 0 0 0 1.81123e+07 1600 1 -1.96764e+10 -2.92287e+09 -1.55519e+10 -2.92287e+09 -4.19079e+09 -2.74221e+09 -1.55519e+10 -2.74221e+09 -1.46328e+10 8.94851e-15 1.81123e+07 - - -3.71243e+06 654602 3.16316e+06 -0.708369 -0.124905 -0.603561 4.56608e+10 1600 -1.94027e-17 4.5615e+10 -1.26129e+08 -5.85141e+06 -1.26129e+08 4.55822e+10 -1.08378e+08 -5.85141e+06 -1.08378e+08 4.57293e+10 9.17168e-14 -4.52843e+06 -3.26466e+06 1.88485e+06 3.16316e+06 -0.659004 -0.297874 -0.603859 4.56587e+10 1600 -1.34824e-18 4.54752e+10 3.89659e+07 -7.38383e+07 3.89659e+07 4.57171e+10 7.66056e+07 -7.38383e+07 7.66056e+07 4.57286e+10 9.127e-14 -6.62725e+06 - - -4.80632e+06 847485 4.0952e+06 0 0 0 1.81123e+07 1600 1 -1.96764e+10 -2.92287e+09 -1.55519e+10 -2.92287e+09 -4.19079e+09 -2.74221e+09 -1.55519e+10 -2.74221e+09 -1.46328e+10 8.94851e-15 1.81123e+07 -4.22661e+06 2.44023e+06 4.0952e+06 -0.0644111 0.111563 3.93309e-17 1.73775e+07 1600 1 -1.68018e+10 -5.75744e+09 -1.36291e+10 -5.75744e+09 -9.28812e+09 -8.51327e+09 -1.36291e+10 -8.51327e+09 -1.41625e+10 8.74508e-15 1.73775e+07 - - -3.26466e+06 1.88485e+06 3.16316e+06 -0.659004 -0.297874 -0.603859 4.56587e+10 1600 -1.34824e-18 4.54752e+10 3.89659e+07 -7.38383e+07 3.89659e+07 4.57171e+10 7.66056e+07 -7.38383e+07 7.66056e+07 4.57286e+10 9.127e-14 -6.62725e+06 -4.22661e+06 2.44023e+06 4.0952e+06 -0.0644111 0.111563 3.93309e-17 1.73775e+07 1600 1 -1.68018e+10 -5.75744e+09 -1.36291e+10 -5.75744e+09 -9.28812e+09 -8.51327e+09 -1.36291e+10 -8.51327e+09 -1.41625e+10 8.74508e-15 1.73775e+07 - - -2.82469e+06 1.63084e+06 1.18715e+06 0.19034 0.0376886 -0.504666 9.13093e+10 1600 0 9.13174e+10 -4.73794e+07 1.8135e+08 -4.73794e+07 9.12319e+10 3.01858e+08 1.8135e+08 3.01858e+08 9.05849e+10 2.67694e-13 -6.9817e+07 -4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56578e+10 -3.05214e+07 1.47808e+06 -3.05214e+07 4.56495e+10 -4.8541e+06 1.47808e+06 -4.8541e+06 4.55995e+10 2.20884e-14 4.8163e+06 - - -2.82469e+06 1.63084e+06 1.18715e+06 0.19034 0.0376886 -0.504666 9.13093e+10 1600 0 9.13174e+10 -4.73794e+07 1.8135e+08 -4.73794e+07 9.12319e+10 3.01858e+08 1.8135e+08 3.01858e+08 9.05849e+10 2.67694e-13 -6.9817e+07 -2.09656e+06 2.49859e+06 1.18715e+06 0.151969 0.147411 -0.578639 9.13366e+10 1600 0 9.11902e+10 -2.31881e+07 3.58683e+08 -2.31881e+07 9.13885e+10 1.07386e+08 3.58683e+08 1.07386e+08 9.05995e+10 2.78025e-13 -4.25455e+07 - - -2.82469e+06 1.63084e+06 1.18715e+06 0.19034 0.0376886 -0.504666 9.13093e+10 1600 0 9.13174e+10 -4.73794e+07 1.8135e+08 -4.73794e+07 9.12319e+10 3.01858e+08 1.8135e+08 3.01858e+08 9.05849e+10 2.67694e-13 -6.9817e+07 -2.30271e+06 1.32947e+06 2.23112e+06 -0.045879 0.0794648 0 9.15743e+10 1600 0 9.17795e+10 1.47943e+08 -2.95352e+08 1.47943e+08 9.16213e+10 -1.71012e+08 -2.95352e+08 -1.71012e+08 9.21342e+10 2.27683e-13 1.95234e+08 - - -4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56578e+10 -3.05214e+07 1.47808e+06 -3.05214e+07 4.56495e+10 -4.8541e+06 1.47808e+06 -4.8541e+06 4.55995e+10 2.20884e-14 4.8163e+06 -2.9724e+06 3.54236e+06 1.68308e+06 -0.120017 -0.203298 -0.283249 4.56924e+10 1600 0 4.56808e+10 -7.47783e+06 -5.79055e+06 -7.47783e+06 4.56548e+10 183550 -5.79055e+06 183550 4.56234e+10 1.51229e-14 2.71349e+07 - - -4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56578e+10 -3.05214e+07 1.47808e+06 -3.05214e+07 4.56495e+10 -4.8541e+06 1.47808e+06 -4.8541e+06 4.55995e+10 2.20884e-14 4.8163e+06 -3.26466e+06 1.88485e+06 3.16316e+06 -0.659004 -0.297874 -0.603859 4.56587e+10 1600 -1.34824e-18 4.57179e+10 -7.08523e+07 5.76269e+07 -7.08523e+07 4.55179e+10 -8.88035e+07 5.76269e+07 -8.88035e+07 4.57547e+10 9.00879e-14 -6.62725e+06 - - -2.09656e+06 2.49859e+06 1.18715e+06 0.151969 0.147411 -0.578639 9.13366e+10 1600 0 9.11902e+10 -2.31881e+07 3.58683e+08 -2.31881e+07 9.13885e+10 1.07386e+08 3.58683e+08 1.07386e+08 9.05995e+10 2.78025e-13 -4.25455e+07 -2.9724e+06 3.54236e+06 1.68308e+06 -0.120017 -0.203298 -0.283249 4.56924e+10 1600 0 4.56808e+10 -7.47783e+06 -5.79055e+06 -7.47783e+06 4.56548e+10 183550 -5.79055e+06 183550 4.56234e+10 1.51229e-14 2.71349e+07 - - -2.09656e+06 2.49859e+06 1.18715e+06 0.151969 0.147411 -0.578639 9.13366e+10 1600 0 9.11902e+10 -2.31881e+07 3.58683e+08 -2.31881e+07 9.13885e+10 1.07386e+08 3.58683e+08 1.07386e+08 9.05995e+10 2.78025e-13 -4.25455e+07 -1.70913e+06 2.03687e+06 2.23112e+06 -0.11234 0.0942646 0 9.15041e+10 1600 0 9.16788e+10 1.36191e+08 -2.31722e+08 1.36191e+08 9.16625e+10 -2.50142e+08 -2.31722e+08 -2.50142e+08 9.2071e+10 2.16969e-13 1.24956e+08 - - -2.9724e+06 3.54236e+06 1.68308e+06 -0.120017 -0.203298 -0.283249 4.56924e+10 1600 0 4.56808e+10 -7.47783e+06 -5.79055e+06 -7.47783e+06 4.56548e+10 183550 -5.79055e+06 183550 4.56234e+10 1.51229e-14 2.71349e+07 -2.42312e+06 2.88776e+06 3.16316e+06 -0.690435 -0.372828 -0.612044 4.56604e+10 1600 0 4.54595e+10 -2.64945e+07 -6.95638e+07 -2.64945e+07 4.57508e+10 6.4528e+07 -6.95638e+07 6.4528e+07 4.575e+10 9.73365e-14 -4.90659e+06 - - -2.30271e+06 1.32947e+06 2.23112e+06 -0.045879 0.0794648 0 9.15743e+10 1600 0 9.17795e+10 1.47943e+08 -2.95352e+08 1.47943e+08 9.16213e+10 -1.71012e+08 -2.95352e+08 -1.71012e+08 9.21342e+10 2.27683e-13 1.95234e+08 -3.26466e+06 1.88485e+06 3.16316e+06 -0.659004 -0.297874 -0.603859 4.56587e+10 1600 -1.34824e-18 4.57179e+10 -7.08523e+07 5.76269e+07 -7.08523e+07 4.55179e+10 -8.88035e+07 5.76269e+07 -8.88035e+07 4.57547e+10 9.00879e-14 -6.62725e+06 - - -2.30271e+06 1.32947e+06 2.23112e+06 -0.045879 0.0794648 0 9.15743e+10 1600 0 9.17795e+10 1.47943e+08 -2.95352e+08 1.47943e+08 9.16213e+10 -1.71012e+08 -2.95352e+08 -1.71012e+08 9.21342e+10 2.27683e-13 1.95234e+08 -1.70913e+06 2.03687e+06 2.23112e+06 -0.11234 0.0942646 0 9.15041e+10 1600 0 9.16788e+10 1.36191e+08 -2.31722e+08 1.36191e+08 9.16625e+10 -2.50142e+08 -2.31722e+08 -2.50142e+08 9.2071e+10 2.16969e-13 1.24956e+08 - - -3.26466e+06 1.88485e+06 3.16316e+06 -0.659004 -0.297874 -0.603859 4.56587e+10 1600 -1.34824e-18 4.57179e+10 -7.08523e+07 5.76269e+07 -7.08523e+07 4.55179e+10 -8.88035e+07 5.76269e+07 -8.88035e+07 4.57547e+10 9.00879e-14 -6.62725e+06 -2.42312e+06 2.88776e+06 3.16316e+06 -0.690435 -0.372828 -0.612044 4.56604e+10 1600 0 4.54595e+10 -2.64945e+07 -6.95638e+07 -2.64945e+07 4.57508e+10 6.4528e+07 -6.95638e+07 6.4528e+07 4.575e+10 9.73365e-14 -4.90659e+06 - - -1.70913e+06 2.03687e+06 2.23112e+06 -0.11234 0.0942646 0 9.15041e+10 1600 0 9.16788e+10 1.36191e+08 -2.31722e+08 1.36191e+08 9.16625e+10 -2.50142e+08 -2.31722e+08 -2.50142e+08 9.2071e+10 2.16969e-13 1.24956e+08 -2.42312e+06 2.88776e+06 3.16316e+06 -0.690435 -0.372828 -0.612044 4.56604e+10 1600 0 4.54595e+10 -2.64945e+07 -6.95638e+07 -2.64945e+07 4.57508e+10 6.4528e+07 -6.95638e+07 6.4528e+07 4.575e+10 9.73365e-14 -4.90659e+06 - - -4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56882e+10 -1.40908e+07 -1.68188e+07 -1.40908e+07 4.56583e+10 -1.59078e+07 -1.68188e+07 -1.59078e+07 4.55787e+10 3.13788e-14 4.8163e+06 -5.18471e+06 2.99339e+06 2.17901e+06 -0.0171684 0.0586576 -0.0397298 -1.83017e+06 1600 1 2.85262e+09 3.43257e+09 3.39323e+09 3.43257e+09 2.24543e+07 -2.1608e+08 3.39323e+09 -2.1608e+08 -5.94477e+08 2.58467e-15 -1.83017e+06 - - -4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56882e+10 -1.40908e+07 -1.68188e+07 -1.40908e+07 4.56583e+10 -1.59078e+07 -1.68188e+07 -1.59078e+07 4.55787e+10 3.13788e-14 4.8163e+06 -2.9724e+06 3.54236e+06 1.68308e+06 -0.120017 -0.203298 -0.283249 4.56924e+10 1600 0 4.56741e+10 1.9163e+06 -2.63568e+07 1.9163e+06 4.56867e+10 -1.44764e+07 -2.63568e+07 -1.44764e+07 4.56023e+10 2.72956e-14 2.71349e+07 - - -4.0047e+06 2.31211e+06 1.68308e+06 -0.118272 -0.12018 -0.23084 4.56701e+10 1600 -6.00786e-18 4.56882e+10 -1.40908e+07 -1.68188e+07 -1.40908e+07 4.56583e+10 -1.59078e+07 -1.68188e+07 -1.59078e+07 4.55787e+10 3.13788e-14 4.8163e+06 -3.26466e+06 1.88485e+06 3.16316e+06 -0.659004 -0.297874 -0.603859 4.56587e+10 1600 -1.34824e-18 4.56915e+10 -8.7796e+07 3.13824e+07 -8.7796e+07 4.55071e+10 -1.05642e+08 3.13824e+07 -1.05642e+08 4.57286e+10 9.20921e-14 -6.62725e+06 - - -5.18471e+06 2.99339e+06 2.17901e+06 -0.0171684 0.0586576 -0.0397298 -1.83017e+06 1600 1 2.85262e+09 3.43257e+09 3.39323e+09 3.43257e+09 2.24543e+07 -2.1608e+08 3.39323e+09 -2.1608e+08 -5.94477e+08 2.58467e-15 -1.83017e+06 -3.84823e+06 4.58614e+06 2.17901e+06 -0.0834525 0.113859 -0.0922569 -1.96489e+07 1600 1 4.76835e+08 5.04687e+09 1.53814e+09 5.04687e+09 4.25411e+09 2.76449e+09 1.53814e+09 2.76449e+09 1.39318e+09 3.13692e-15 -1.96489e+07 - - -5.18471e+06 2.99339e+06 2.17901e+06 -0.0171684 0.0586576 -0.0397298 -1.83017e+06 1600 1 2.85262e+09 3.43257e+09 3.39323e+09 3.43257e+09 2.24543e+07 -2.1608e+08 3.39323e+09 -2.1608e+08 -5.94477e+08 2.58467e-15 -1.83017e+06 -4.22661e+06 2.44023e+06 4.0952e+06 -0.0644111 0.111563 3.93309e-17 1.73775e+07 1600 1 -1.6204e+10 -6.51291e+09 -1.36291e+10 -6.51291e+09 -8.46455e+09 -8.51325e+09 -1.36291e+10 -8.51325e+09 -1.41625e+10 8.89851e-15 1.73775e+07 - - -2.9724e+06 3.54236e+06 1.68308e+06 -0.120017 -0.203298 -0.283249 4.56924e+10 1600 0 4.56741e+10 1.9163e+06 -2.63568e+07 1.9163e+06 4.56867e+10 -1.44764e+07 -2.63568e+07 -1.44764e+07 4.56023e+10 2.72956e-14 2.71349e+07 -3.84823e+06 4.58614e+06 2.17901e+06 -0.0834525 0.113859 -0.0922569 -1.96489e+07 1600 1 4.76835e+08 5.04687e+09 1.53814e+09 5.04687e+09 4.25411e+09 2.76449e+09 1.53814e+09 2.76449e+09 1.39318e+09 3.13692e-15 -1.96489e+07 - - -2.9724e+06 3.54236e+06 1.68308e+06 -0.120017 -0.203298 -0.283249 4.56924e+10 1600 0 4.56741e+10 1.9163e+06 -2.63568e+07 1.9163e+06 4.56867e+10 -1.44764e+07 -2.63568e+07 -1.44764e+07 4.56023e+10 2.72956e-14 2.71349e+07 -2.42312e+06 2.88776e+06 3.16316e+06 -0.690435 -0.372828 -0.612044 4.56604e+10 1600 0 4.54514e+10 -4.2775e+07 -8.43876e+07 -4.2775e+07 4.57234e+10 3.81818e+07 -8.43876e+07 3.81818e+07 4.5725e+10 9.38237e-14 -4.90659e+06 - - -3.84823e+06 4.58614e+06 2.17901e+06 -0.0834525 0.113859 -0.0922569 -1.96489e+07 1600 1 4.76835e+08 5.04687e+09 1.53814e+09 5.04687e+09 4.25411e+09 2.76449e+09 1.53814e+09 2.76449e+09 1.39318e+09 3.13692e-15 -1.96489e+07 -3.13711e+06 3.73866e+06 4.0952e+06 -0.291372 0.244491 2.21699e-16 -1.06766e+07 1600 1 -1.77671e+10 -3.19649e+09 -1.66819e+09 -3.19649e+09 -1.75129e+10 -1.52152e+10 -1.66819e+09 -1.52152e+10 -1.68353e+10 7.82202e-15 -1.06766e+07 - - -3.26466e+06 1.88485e+06 3.16316e+06 -0.659004 -0.297874 -0.603859 4.56587e+10 1600 -1.34824e-18 4.56915e+10 -8.7796e+07 3.13824e+07 -8.7796e+07 4.55071e+10 -1.05642e+08 3.13824e+07 -1.05642e+08 4.57286e+10 9.20921e-14 -6.62725e+06 -4.22661e+06 2.44023e+06 4.0952e+06 -0.0644111 0.111563 3.93309e-17 1.73775e+07 1600 1 -1.6204e+10 -6.51291e+09 -1.36291e+10 -6.51291e+09 -8.46455e+09 -8.51325e+09 -1.36291e+10 -8.51325e+09 -1.41625e+10 8.89851e-15 1.73775e+07 - - -3.26466e+06 1.88485e+06 3.16316e+06 -0.659004 -0.297874 -0.603859 4.56587e+10 1600 -1.34824e-18 4.56915e+10 -8.7796e+07 3.13824e+07 -8.7796e+07 4.55071e+10 -1.05642e+08 3.13824e+07 -1.05642e+08 4.57286e+10 9.20921e-14 -6.62725e+06 -2.42312e+06 2.88776e+06 3.16316e+06 -0.690435 -0.372828 -0.612044 4.56604e+10 1600 0 4.54514e+10 -4.2775e+07 -8.43876e+07 -4.2775e+07 4.57234e+10 3.81818e+07 -8.43876e+07 3.81818e+07 4.5725e+10 9.38237e-14 -4.90659e+06 - - -4.22661e+06 2.44023e+06 4.0952e+06 -0.0644111 0.111563 3.93309e-17 1.73775e+07 1600 1 -1.6204e+10 -6.51291e+09 -1.36291e+10 -6.51291e+09 -8.46455e+09 -8.51325e+09 -1.36291e+10 -8.51325e+09 -1.41625e+10 8.89851e-15 1.73775e+07 -3.13711e+06 3.73866e+06 4.0952e+06 -0.291372 0.244491 2.21699e-16 -1.06766e+07 1600 1 -1.77671e+10 -3.19649e+09 -1.66819e+09 -3.19649e+09 -1.75129e+10 -1.52152e+10 -1.66819e+09 -1.52152e+10 -1.68353e+10 7.82202e-15 -1.06766e+07 - - -2.42312e+06 2.88776e+06 3.16316e+06 -0.690435 -0.372828 -0.612044 4.56604e+10 1600 0 4.54514e+10 -4.2775e+07 -8.43876e+07 -4.2775e+07 4.57234e+10 3.81818e+07 -8.43876e+07 3.81818e+07 4.5725e+10 9.38237e-14 -4.90659e+06 -3.13711e+06 3.73866e+06 4.0952e+06 -0.291372 0.244491 2.21699e-16 -1.06766e+07 1600 1 -1.77671e+10 -3.19649e+09 -1.66819e+09 -3.19649e+09 -1.75129e+10 -1.52152e+10 -1.66819e+09 -1.52152e+10 -1.68353e+10 7.82202e-15 -1.06766e+07 +3.15254e+06 3.75706e+06 -0.208402 -0.424846 4.7355e+10 1600 0 4.73355e+10 -599213 -599213 4.73009e+10 8.64366e-15 1.1692e+09 +4.07397e+06 4.85517e+06 0.247802 -0.199201 2.26254e+09 1600 1 -1.58507e+10 -3.16023e+10 -3.16023e+10 6.66552e+10 2.59832e-14 1.23724e+09 diff --git a/tests/chunk_lithostatic_pressure_2d_initial_topography/statistics b/tests/chunk_lithostatic_pressure_2d_initial_topography/statistics index 7c4ebdcbbd0..162d6c2a140 100644 --- a/tests/chunk_lithostatic_pressure_2d_initial_topography/statistics +++ b/tests/chunk_lithostatic_pressure_2d_initial_topography/statistics @@ -16,4 +16,4 @@ # 16: Average pressure (Pa) # 17: Maximal pressure (Pa) # 18: Visualization file name -0 0.000000000000e+00 0.000000000000e+00 8 402 125 125 0 11 21 23 22 3.09329044e-01 1.86122381e+00 -1.96488886e+07 3.29000931e+10 9.15743340e+10 output-chunk_lithostatic_pressure_2d/solution/solution-00000 +0 0.000000000000e+00 0.000000000000e+00 4 59 25 25 0 8 18 20 38 2.74528165e-01 4.28340212e-01 -2.13978640e+09 3.92771590e+10 9.24905403e+10 output-chunk_lithostatic_pressure_2d_initial_topography/solution/solution-00000 From 1c865b534c9309ffd33925ed2f4ddf6c345a24f7 Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Thu, 14 Dec 2023 11:57:53 +0100 Subject: [PATCH 17/22] Add change log --- doc/modules/changes/20231214_glerum | 5 +++++ 1 file changed, 5 insertions(+) create mode 100644 doc/modules/changes/20231214_glerum diff --git a/doc/modules/changes/20231214_glerum b/doc/modules/changes/20231214_glerum new file mode 100644 index 00000000000..b72905b8842 --- /dev/null +++ b/doc/modules/changes/20231214_glerum @@ -0,0 +1,5 @@ + Improved: The geometry function point_is_in_domain + now also works for meshes that include initial topography + and/or mesh deformation for box and chunk geometries. +
+ (Anne Glerum, 2023/12/14) From 1fbed4f7385208562b05ce1b0ab291d32e9928cf Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Thu, 14 Dec 2023 14:24:18 +0100 Subject: [PATCH 18/22] Add utilities function --- include/aspect/utilities.h | 10 ++++++++++ source/utilities.cc | 31 +++++++++++++++++++++++++++++-- 2 files changed, 39 insertions(+), 2 deletions(-) diff --git a/include/aspect/utilities.h b/include/aspect/utilities.h index 13578eb64c6..85d27c4296e 100644 --- a/include/aspect/utilities.h +++ b/include/aspect/utilities.h @@ -56,6 +56,16 @@ namespace aspect using namespace dealii; using namespace dealii::Utilities; + /** Given an + * + */ + template + bool + point_is_in_triangulation(const Mapping &mapping, + const parallel::distributed::Triangulation &triangulation, + const Point point, + const MPI_Comm mpi_communicator); + /** * Given an array @p values, consider three cases: * - If it has size @p N, return the original array. diff --git a/source/utilities.cc b/source/utilities.cc index 480f1d7e23b..cacec8c9f01 100644 --- a/source/utilities.cc +++ b/source/utilities.cc @@ -39,8 +39,7 @@ #include #include #include - - +#include #include #include @@ -214,6 +213,34 @@ namespace aspect + template + bool + point_is_in_triangulation(const Mapping &mapping, + const parallel::distributed::Triangulation &triangulation, + const Point point, + const MPI_Comm mpi_communicator) + { + // Try to find the cell around the given point. + bool cell_found = false; + std::pair::active_cell_iterator, + Point> it = + GridTools::find_active_cell_around_point<>(mapping, triangulation, point); + + // If we found the correct cell on this MPI process, we have found the right cell. + if (it.first.state() == IteratorState::valid && it.first->is_locally_owned()) + cell_found = true; + + // Compute how many processes found the cell. + const int n_procs_cell_found = Utilities::MPI::sum(cell_found ? 1 : 0, mpi_communicator); + + // If at least one process found the cell, the point is in the triangulation. + if (n_procs_cell_found > 0) + return true; + else + return false; + } + + template Table<2,T> parse_input_table (const std::string &input_string, From 8e15a1139ea2c3b35ff2eeb55178f2e3c524ef3e Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Thu, 14 Dec 2023 15:57:06 +0100 Subject: [PATCH 19/22] Space test output --- tests/airy_isostasy_initial_topo.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/airy_isostasy_initial_topo.cc b/tests/airy_isostasy_initial_topo.cc index 3f194ab00e9..6576120d3b1 100644 --- a/tests/airy_isostasy_initial_topo.cc +++ b/tests/airy_isostasy_initial_topo.cc @@ -53,8 +53,8 @@ namespace aspect const Point point_not_in_domain_2(2.0, 0.01); std::ostringstream screen_text; - screen_text << std::boolalpha << this->get_geometry_model().point_is_in_domain(point_in_domain) ; - screen_text << std::boolalpha << this->get_geometry_model().point_is_in_domain(point_not_in_domain); + screen_text << std::boolalpha << this->get_geometry_model().point_is_in_domain(point_in_domain) << " "; + screen_text << std::boolalpha << this->get_geometry_model().point_is_in_domain(point_not_in_domain) << " "; screen_text << std::boolalpha << this->get_geometry_model().point_is_in_domain(point_not_in_domain_2); return std::pair ("Points lie in domain: ", From 2ce8227932f863e3e8f26dabfb750e8e631fe73f Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Thu, 14 Dec 2023 18:06:18 +0100 Subject: [PATCH 20/22] Update utilities function --- include/aspect/utilities.h | 22 ++++++++------ source/utilities.cc | 59 ++++++++++++++++++++++---------------- 2 files changed, 47 insertions(+), 34 deletions(-) diff --git a/include/aspect/utilities.h b/include/aspect/utilities.h index 85d27c4296e..aafaf1d3e23 100644 --- a/include/aspect/utilities.h +++ b/include/aspect/utilities.h @@ -56,15 +56,6 @@ namespace aspect using namespace dealii; using namespace dealii::Utilities; - /** Given an - * - */ - template - bool - point_is_in_triangulation(const Mapping &mapping, - const parallel::distributed::Triangulation &triangulation, - const Point point, - const MPI_Comm mpi_communicator); /** * Given an array @p values, consider three cases: @@ -366,6 +357,19 @@ namespace aspect + /** + * Given a point @p point, find out if any of the MPI + * processes own the cell in which this point lies. If + * not, the point lies outside the @p triangulation. + */ + template + bool + point_is_in_triangulation(const Mapping &mapping, + const parallel::distributed::Triangulation &triangulation, + const Point &point, + const MPI_Comm mpi_communicator); + + namespace Coordinates { diff --git a/source/utilities.cc b/source/utilities.cc index cacec8c9f01..6665c5e22ac 100644 --- a/source/utilities.cc +++ b/source/utilities.cc @@ -213,32 +213,7 @@ namespace aspect - template - bool - point_is_in_triangulation(const Mapping &mapping, - const parallel::distributed::Triangulation &triangulation, - const Point point, - const MPI_Comm mpi_communicator) - { - // Try to find the cell around the given point. - bool cell_found = false; - std::pair::active_cell_iterator, - Point> it = - GridTools::find_active_cell_around_point<>(mapping, triangulation, point); - // If we found the correct cell on this MPI process, we have found the right cell. - if (it.first.state() == IteratorState::valid && it.first->is_locally_owned()) - cell_found = true; - - // Compute how many processes found the cell. - const int n_procs_cell_found = Utilities::MPI::sum(cell_found ? 1 : 0, mpi_communicator); - - // If at least one process found the cell, the point is in the triangulation. - if (n_procs_cell_found > 0) - return true; - else - return false; - } template @@ -756,6 +731,34 @@ namespace aspect + template + bool + point_is_in_triangulation(const Mapping &mapping, + const parallel::distributed::Triangulation &triangulation, + const Point &point, + const MPI_Comm mpi_communicator) + { + // Try to find the cell around the given point. + bool cell_found = false; + std::pair::active_cell_iterator, + Point> it = + GridTools::find_active_cell_around_point<>(mapping, triangulation, point); + + // If we found the correct cell on this MPI process, we have found the right cell. + if (it.first.state() == IteratorState::valid && it.first->is_locally_owned()) + cell_found = true; + + // Compute how many processes found the cell. + const int n_procs_cell_found = Utilities::MPI::sum(cell_found ? 1 : 0, mpi_communicator); + // If at least one process found the cell, the point is in the triangulation. + if (n_procs_cell_found > 0) + return true; + else + return false; + } + + + namespace Coordinates { @@ -3455,6 +3458,12 @@ namespace aspect const dealii::Point<2> &point); \ \ template \ + bool point_is_in_triangulation(const Mapping &mapping, \ + const parallel::distributed::Triangulation &triangulation, \ + const Point &point, \ + const MPI_Comm mpi_communicator); \ + \ + template \ double signed_distance_to_polygon(const std::vector> &pointList, \ const dealii::Point<2> &point); \ \ From f923ae1adbf1505c716b602bd932e7c15fa5a561 Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Thu, 14 Dec 2023 18:12:21 +0100 Subject: [PATCH 21/22] Use function and simplify comment --- source/geometry_model/box.cc | 53 ++++++++++------------------------ source/geometry_model/chunk.cc | 36 +++++++---------------- 2 files changed, 26 insertions(+), 63 deletions(-) diff --git a/source/geometry_model/box.cc b/source/geometry_model/box.cc index 6d17441884a..4d796ee6b11 100644 --- a/source/geometry_model/box.cc +++ b/source/geometry_model/box.cc @@ -23,6 +23,7 @@ #include #include #include +#include #include #include @@ -317,51 +318,27 @@ namespace aspect bool Box::point_is_in_domain(const Point &point) const { - // If mesh deformation occurs over time, and we're past the - // first timestep (during which no mesh deformation is applied), - // we can no longer use the original extents of the model domain - // to check whether the point lies in the domain. - // Initial mesh deformation is applied during the first - // timestep, but this function might be called before the - // simulator is initialized and thus before we can ask whether initial - // or time-dependent mesh deformation is applied. - // Therefore, when any mesh deformation is enabled, - // we loop over all the current grid cells to see if the point - // lies within the domain. - // Do note that this function can be called before mesh deformation - // is applied in the first timestep, and therefore there is no guarantee - // that the point will lie in the domain after initial mesh deformation. - // For example, boundary traction plugins are initialized before - // mesh deformation plugins and the initial_lithostatic_pressure - // plugin calls this function while the simulator is being initialized. - // However, if the simulator has not been initialized, - // no mesh deformation has been applied and we can use the geometry model's - // extents to check whether a point lies in the domain. + // If mesh deformation is enabled, we have to loop over all the current + // grid cells to see if the given point lies in the domain. + // If mesh deformation is not enabled, or has not happened yet, + // we can use the global extents of the model domain with or without + // initial topography instead. + // This function only checks if the given point lies in the domain + // in its current shape at the current time. It can be called before + // mesh deformation is applied in the first timestep (e.g., by the boundary + // traction plugins), and therefore there is no guarantee + // that the point will still lie in the domain after initial mesh deformation. if (this->get_parameters().mesh_deformation_enabled && this->simulator_is_past_initialization()) { - bool cell_found = false; - std::pair::active_cell_iterator, - Point> it = - GridTools::find_active_cell_around_point<> (this->get_mapping(), this->get_triangulation(), point); - - // If the cell the point is in is on this processor, we have found the right cell. - if (it.first.state() == IteratorState::valid && it.first->is_locally_owned()) - cell_found = true; - - // Compute how many processes found the cell. - const int n_procs_cell_found = Utilities::MPI::sum(cell_found ? 1 : 0, this->get_mpi_communicator()); - - // If at least one process found the cell, the point is in the domain. - if (n_procs_cell_found > 0) - return true; - else - return false; + return Utilities::point_is_in_triangulation(this->get_mapping(), + this->get_triangulation(), + point, + this->get_mpi_communicator()); } // Without mesh deformation enabled, it is much cheaper to check whether the point lies in the domain. else { - // The maximal extents of the unperturbed box domain. Point max_point = extents+box_origin; diff --git a/source/geometry_model/chunk.cc b/source/geometry_model/chunk.cc index 36049380d0c..d4195f4fb5c 100644 --- a/source/geometry_model/chunk.cc +++ b/source/geometry_model/chunk.cc @@ -697,34 +697,20 @@ namespace aspect bool Chunk::point_is_in_domain(const Point &point) const { - // If mesh deformation has possibly been applied, - // e.g. initial mesh deformation after the simulator - // has been initialized or time-dependent mesh deformation - // after timestep 0, we cannot use information from the - // geometry model to determine whether the point lies in - // the model domain. Instead, we loop over all cells - // to check whether the point lies in one of them. + // If mesh deformation is enabled, we have to loop over all the current + // grid cells to see if the given point lies in the domain. + // If mesh deformation is not enabled, or has not happened yet, + // we can use the global extents of the model domain with or without + // initial topography instead. + // This function only checks if the given point lies in the domain + // in its current shape at the current time. It can be called before + // mesh deformation is applied in the first timestep (e.g., by the boundary + // traction plugins), and therefore there is no guarantee + // that the point will still lie in the domain after initial mesh deformation. if (this->get_parameters().mesh_deformation_enabled && this->simulator_is_past_initialization()) { - bool cell_found = false; - std::pair::active_cell_iterator, - Point> - it = - GridTools::find_active_cell_around_point<>(this->get_mapping(), this->get_triangulation(), point); - - // If the cell the point is in is on this processor, we have found the right cell. - if (it.first.state() == IteratorState::valid && it.first->is_locally_owned()) - cell_found = true; - - // Compute how many processes found the cell. - const int n_procs_cell_found = Utilities::MPI::sum(cell_found ? 1 : 0, this->get_mpi_communicator()); - - // If at least one process found the cell, the point is in the domain. - if (n_procs_cell_found > 0) - return true; - else - return false; + return Utilities::point_is_in_triangulation(this->get_mapping(), this->get_triangulation(), point, this->get_mpi_communicator()); } // Without mesh deformation enabled, it is much cheaper to check whether the point lies in the domain. else From 4b997e0c94d7b391ae4a31b7caa59d7bea35cef1 Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Thu, 14 Dec 2023 21:02:51 +0100 Subject: [PATCH 22/22] Update test output --- tests/airy_isostasy_initial_topo/screen-output | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/airy_isostasy_initial_topo/screen-output b/tests/airy_isostasy_initial_topo/screen-output index e13824c064d..7211f5ab528 100644 --- a/tests/airy_isostasy_initial_topo/screen-output +++ b/tests/airy_isostasy_initial_topo/screen-output @@ -43,7 +43,7 @@ Number of mesh deformation degrees of freedom: 4,162 RMS, max velocity: 1.76e+06 m/year, 2.3e+06 m/year Pressure min/avg/max: -1.132 Pa, 5.012 Pa, 10 Pa Topography min/max: 0 m, 0.05 m - Points lie in domain: truefalsefalse + Points lie in domain: true false false *** Timestep 1: t=1.69957e-09 years, dt=1.69957e-09 years Solving mesh displacement system... 8 iterations. @@ -56,7 +56,7 @@ Number of mesh deformation degrees of freedom: 4,162 RMS, max velocity: 1.62e+06 m/year, 2.11e+06 m/year Pressure min/avg/max: -1.119 Pa, 5.004 Pa, 10.01 Pa Topography min/max: -4.848e-05 m, 0.04609 m - Points lie in domain: truefalsefalse + Points lie in domain: true false false *** Timestep 2: t=3.5e-09 years, dt=1.80043e-09 years Solving mesh displacement system... 8 iterations. @@ -69,7 +69,7 @@ Number of mesh deformation degrees of freedom: 4,162 RMS, max velocity: 1.53e+06 m/year, 1.99e+06 m/year Pressure min/avg/max: -0.8445 Pa, 5.003 Pa, 10.02 Pa Topography min/max: -9.971e-05 m, 0.04229 m - Points lie in domain: truefalsefalse + Points lie in domain: true false false Termination requested by criterion: end time