Chunk initial topography #1705
Chunk initial topography #1705
Conversation
There was a problem hiding this comment.
The current version does not run, but instead creates the following error:
--------------------------------------------------------
An error occurred in line <6280> of file </home/rengas/Software/dealii/include/deal.II/numerics/vector_tools.templates.h> in function
void dealii::VectorTools::compute_nonzero_normal_flux_constraints(const DoFHandlerType<dim, spacedim>&, unsigned int, const std::set<unsigned char>&, typename dealii::FunctionMap<spacedim>::type&, dealii::ConstraintMatrix&, const dealii::Mapping<dim, spacedim>&) [with int dim = 2; int spacedim = 2; DoFHandlerType = dealii::DoFHandler; typename dealii::FunctionMap<spacedim>::type = std::map<unsigned char, const dealii::Function<2>*, std::less<unsigned char>, std::allocator<std::pair<const unsigned char, const dealii::Function<2>*> > >]
The violated condition was:
std::fabs(determinant (t)) > 1e-3
Additional information:
Found a set of normal vectors that are nearly collinear.
Stacktrace:
-----------
#0 /home/rengas/Software/deal.II-dev/lib/libdeal_II.g.so.9.0.0-pre: void dealii::VectorTools::compute_nonzero_normal_flux_constraints<2, 2, dealii::DoFHandler>(dealii::DoFHandler<2, 2> const&, unsigned int, std::set<unsigned char, std::less<unsigned char>, std::allocator<unsigned char> > const&, dealii::FunctionMap<2, double>::type&, dealii::ConstraintMatrix&, dealii::Mapping<2, 2> const&)
#1 /home/rengas/Software/deal.II-dev/lib/libdeal_II.g.so.9.0.0-pre: void dealii::VectorTools::compute_no_normal_flux_constraints<2, 2, dealii::DoFHandler>(dealii::DoFHandler<2, 2> const&, unsigned int, std::set<unsigned char, std::less<unsigned char>, std::allocator<unsigned char> > const&, dealii::ConstraintMatrix&, dealii::Mapping<2, 2> const&)
#2 aspect: aspect::Simulator<2>::setup_dofs()
#3 aspect: aspect::Simulator<2>::run()
#4 aspect: main
--------------------------------------------------------
I was hoping to find a clue about the problem by reviewing the code, but it is not obvious to me where the problem is. One clue is that if I remove the tangential velocity boundaries the model runs, and it does not seem to produce any topography independent on the size of topography I prescribe in the data file. Need to take a closer look in the next days, but feel free to go ahead and fix my comments, maybe you notice something that is not right.
| @@ -134,6 +132,11 @@ namespace aspect | |||
| double depth(const Point<dim> &position) const; | |||
|
|
|||
| virtual | |||
| double depth_wrt_topo(const Point<dim> &position) const; | |||
| @@ -269,6 +272,13 @@ namespace aspect | |||
| public: | |||
| ChunkGeometry(); | |||
|
|
|||
| /** | |||
| * An initialization function necessary for to make sure that the | |||
|
|
||
| virtual | ||
| Point<dim> | ||
| push_forward_topo(const Point<dim> &chart_point) const; |
There was a problem hiding this comment.
Please add documentation to all of the functions above. You can copy a lot from EllipsoidalChunk
| @@ -68,6 +68,12 @@ namespace aspect | |||
| double | |||
| value (const Point<dim-1> &p) const; | |||
|
|
|||
| /** | |||
| * Return the gradient of the surface topography for a given position along | |||
| * along the surface. | |||
| * Return the value of the elevation at the given point. | ||
| */ | ||
| // virtual | ||
| // Tensor<1,dim-1> vector_gradient (const Point<dim> &p) const = 0; |
There was a problem hiding this comment.
I think this function is still here and unnecessary?
source/geometry_model/chunk.cc
Outdated
| Point<dim> topor_phi_theta = r_phi_theta; | ||
| topor_phi_theta[0] = topo_radius; | ||
| return topor_phi_theta; | ||
| // return r_phi_theta; |
source/geometry_model/chunk.cc
Outdated
| // Convert latitude to colatitude | ||
| if (dim == 3) | ||
| surface_point[1] = 0.5*numbers::PI - surface_point[1]; | ||
| topography = topo->value(surface_point); |
There was a problem hiding this comment.
make to const double topography
source/geometry_model/chunk.cc
Outdated
| Point<dim> r_phi_theta = topor_phi_theta; | ||
| r_phi_theta[0] = radius; | ||
| return r_phi_theta; | ||
| // return topor_phi_theta; |
source/geometry_model/chunk.cc
Outdated
| @@ -291,6 +542,7 @@ namespace aspect | |||
| } | |||
|
|
|||
| // Attach shell boundary objects to the curved inner and outer boundaries | |||
| // TODO is this boundary still correct for a displaced mesh? | |||
There was a problem hiding this comment.
This is not used with a displaced mesh (only for deal.II >9.0). Or do you mean in case of a free surface unrelated to initial topography?
| // depth is therefore always positive | ||
| const double outer_radius = manifold.get_radius(position); | ||
| const Point<dim> rtopo_phi_theta = manifold.pull_back_sphere(position); | ||
| return std::max(0.0, outer_radius - rtopo_phi_theta[0]); |
There was a problem hiding this comment.
This way you cut off all depth values above the unmodified surface, is that what you want? depth_wrt_topo suggests to me to use the full manifold.pull_back, or do I misunderstand?
There was a problem hiding this comment.
Depth here is computed as the difference between the radius of modified surface (given by get_radius) and the modified point (given by pull_back_sphere).
|
@gassmoeller Thanks for testing! I'll have a look today, now that I can test both versions. Update: Enabling the As we will not use initial topography with the old boundary objects, I’ll leave this as is. The |
|
Great, that sounds like you found the problem 👍. Let us stay consistent and use the signal to work around it. I do not quite remember the logic in |
|
/run-tests |
|
Alright, I'll try and use the signal and address your comments. I ran into the HyperShellBoundary error when I enabled initial topography with deal.II < 9.0 (because I knew this used to work). |
anne-glerum
force-pushed
the
chunk_initial_topo
branch
from
0aa8e8b
to
b1e97f8
Compare
|
@gassmoeller This functionality now also works for dealii 9-pre. Three things:
|
anne-glerum
force-pushed
the
chunk_initial_topo
branch
from
aa089c5
to
6faab39
Compare
|
Hi @gassmoeller @naliboff , I've rebased, updated and simplified this pr and I think it's ready for another review. Repeating the test of which I posted figures above, now results in a flow field from high to low topography (velocity field and vectors): |
|
For reference, I'm posting some snapshots of 1 timestep chunk models with varying resolution/refinement/topography using the updates in this pr. Models are based on the test No topography, only global refinement (bottom) or with adaptive refinement (top), refinement increases to the right: |
|
@anne-glerum - Thanks for updating all of this! I compiled the branch to run a few tests, but unfortunately both the 2D chunk and the models I made are failing at this line in Is the test successfully running on your end? If so, perhaps it is a deal.II version issue? |
|
Hi @naliboff , yes the test works for me. What error message do you get? |
|
@anne-glerum The error message is just a floating point exception, but clearly no division by zero in the failing line. Perhaps the initial data file model time may not be initialized yet? Let me know if digging deeper would be helpful! FYI, I built the branch with deal.II v9.0.1. |
|
@naliboff I've switched the conditional statements. Only in the very specific case of Initial topography specified by an ascii data table, we should get to these lines before time is set to something other than NaN. Now we should hit the condition that checks for NaN first and hopefully that solves the problem. If not, some more info would be helpful and I'll try to reproduce the error. |
anne-glerum
force-pushed
the
chunk_initial_topo
branch
from
e82fc79
to
90e7c10
Compare
|
@anne-glerum - hmmm, still getting the same error with the last commit. I'll dig into the debugger tomorrow and see if I can pinpoint a solution on my end. |
anne-glerum
force-pushed
the
chunk_initial_topo
branch
from
ed574fa
to
e3b93a8
Compare
|
new is_time_initialized parameter could be updated after new function discussed in #2529. Actually, we could use get_timestep_number() == numbers::invalid_unsigned_int |
| * Return the value of the elevation at the given point. | ||
| */ | ||
| // virtual | ||
| // Tensor<1,dim-1> vector_gradient (const Point<dim> &p) const = 0; |
There was a problem hiding this comment.
I think this function is still here and unnecessary?
| private: | ||
| // The minimum longitude of the domain | ||
| double point1_lon; | ||
| double inner_radius; | ||
|
|
||
| double max_depth; |
There was a problem hiding this comment.
Please add comments for max_depth and inner_radius
include/aspect/utilities.h
Outdated
| /** | ||
| * Returns the gradient of the function based on the bilinear | ||
| * interpolation of the data (velocity, temperature, etc. - according | ||
| * to the used plugin) in Cartesian coordinates. TOOD Cartesian?? |
source/geometry_model/chunk.cc
Outdated
| @@ -32,6 +33,10 @@ | |||
| #include <deal.II/grid/tria_boundary_lib.h> | |||
| #endif | |||
|
|
|||
| #if DEAL_II_VERSION_GTE(9,0,0) | |||
There was a problem hiding this comment.
No need for DEAL_II_VERSION_GTE(9,0,0) anymore, because we now require 9.0.
source/geometry_model/chunk.cc
Outdated
|
|
||
| // prior to deal.II 9, we do not apply initial topography | ||
| #if !DEAL_II_VERSION_GTE(9,0,0) | ||
| const double phi = chart_point[1]; // Longitude |
There was a problem hiding this comment.
This case no longer happens, remove the if and the line.
| @@ -1829,7 +1838,6 @@ namespace aspect | |||
| } | |||
|
|
|||
|
|
|||
|
|
|||
| @@ -2165,7 +2173,11 @@ namespace aspect | |||
| const Point<dim> &position, | |||
| const unsigned int component) const | |||
| { | |||
| if (this->get_time() - first_data_file_model_time >= 0.0) | |||
| // For initial ascii data topography, we need access to the data before get_time() is set | |||
| if ( (dynamic_cast<const GeometryModel::Chunk<dim>*>(&this->get_geometry_model()) != nullptr && | |||
There was a problem hiding this comment.
please add a comment why these geometries are special
| return time_weight * gradients + (1 - time_weight) * old_gradients; | ||
| } | ||
| else | ||
| return Tensor<1,dim-1>(); |
There was a problem hiding this comment.
This way you silently return 0.0 for other geometries than the Chunk? Please make this an AssertThrow(false, ...) so that nobody accidentally calls this function for other geometries.
There was a problem hiding this comment.
I think return this Tensor makes sense for the other condition (this->get_time() - first_data_file_model_time >= 0.0). Then if we're not far enough in time yet, 0 is returned. The same is done for get_data_component.
source/geometry_model/chunk.cc
Outdated
|
|
||
| // Grab lon,lat coordinates | ||
| Point<dim-1> surface_point; | ||
| for (unsigned int d=0; d<dim-1; d++) |
source/geometry_model/chunk.cc
Outdated
| @@ -186,7 +340,6 @@ namespace aspect | |||
| } | |||
|
|
|||
|
|
|||
|
|
|||
| #if DEAL_II_VERSION_GTE(9,0,0) | |||
|
@gassmoeller I've addressed the comments. |
|
This PR has been superseded by #3039, so I'm closing it. |
@gassmoeller
Hi Rene,
do you think you could give this geometry model a try? I've added a test that you could run. Also, right now, that test fails on my laptop because I put in an AssertThrow for older deal.II versions that try to add non-zero initial topography. What would be the best way to go around this?
I'd like a more compact way to put in the derivative gradient, but that's for tomorrow :)
Thank you!