Fastscape coupling [WIP]

[Djneu]

This is an initial pull request for review and discussion of the code. In the future, I plan to squash some of the commits and split the pull request into a few different pieces to hopefully make it easier to merge and a bit less messy. @anne-glerum @naliboff

Before your first pull request:

• [ ] I have read the guidelines in our CONTRIBUTING.md document.

For all pull requests:

• [ ] I have followed the instructions for indenting my code.

For new features/models or changes of existing features:

• [ ] I have tested my new feature locally to ensure it is correct.
• [ ] I have created a testcase for the new feature/benchmark in the tests/ directory.
• [ ] I have added a changelog entry in the doc/modules/changes directory that will inform other users of my change.

[Djneu]

Hi all, Just thought I'd add a bit of info since this pull request is much messier than I thought it would be. There are really only ~5 changes that are needed for the pr, the other changes are mostly from a bad indent that reformatted many files for some reason (e.g., all the files in mesh refinement).

1. Fastscape.cc and fastscape.h
2. Changes to box.cc and box.h to add a few functions
3. Cmakelists.txt
4. Named_VTK.f90, which is a file that needs to be moved to the fastscape source directory on installation
5. fastscape_eroding_box.prm, a new cookbook to test the installation.
6. Changes to free_surface and mesh_deformation/interface that relocates the free surface stabilization factor, though this is open in a separate pull request (#4350).

[anne-glerum]

Hi Derek, I know you were thinking of opening up a new PR to get rid of the files that were included by accident here (as you explain in your last comment), but I finished my review here to keep all the comments on the same files in the same place. Also, there might be a solution where we can keep one PR? I can think of checking out the old versions of the specific files and committing them. By squashing at least the edits and their undoing into one commit, maybe the changes will disappear. But maybe that's more hassle than it's worth. @tjhei would you have a solution?

[tjhei]

would you have a solution?

I would do the following:

1. delete all files that are not necessary with git rm, commit
2. run make indent, commit
3. do a git rebase -i and squash all your commits into one
4. force push

[anne-glerum]

For documenting to address in the future: Running a model that has an origin that is not at (0,0) leads to a crash. I have not investigated it further, but it will fail after writing the Initial VTK file in t1. Error message:

*** Error in `/home/projects/bbp00039/fastscape_build_files/aspect-22feb22-dealii-10.0-pre-fastscape-esther': free(): invalid size: 0x0000000005bd4070 ***
======= Backtrace: =========
/lib64/libc.so.6(+0x81329)[0x2aaad167b329]
/home/projects/bbp00039/fastscape_build_files/aspect-22feb22-dealii-10.0-pre-fastscape-esther(_ZNK6aspect15MeshDeformation9FastScapeILi2EE40compute_velocity_constraints_on_boundaryERKN6dealii10DoFHandlerILi2ELi2EEERNS3_17AffineConstraintsIdEERKSt3setIjSt4lessIjESaIjEE+0x2a65)[0x1387045]
/home/projects/bbp00039/fastscape_build_files/aspect-22feb22-dealii-10.0-pre-fastscape-esther(_ZN6aspect15MeshDeformation22MeshDeformationHandlerILi2EE16make_constraintsEv+0x7fd)[0x13a02cd]
/home/projects/bbp00039/fastscape_build_files/aspect-22feb22-dealii-10.0-pre-fastscape-esther(_ZN6aspect15MeshDeformation22MeshDeformationHandlerILi2EE7executeEv+0xc0)[0x13b1ab0]
/home/projects/bbp00039/fastscape_build_files/aspect-22feb22-dealii-10.0-pre-fastscape-esther(_ZN6aspect9SimulatorILi2EE14solve_timestepEv+0x87)[0xa86f57]
/home/projects/bbp00039/fastscape_build_files/aspect-22feb22-dealii-10.0-pre-fastscape-esther(_ZN6aspect9SimulatorILi2EE3runEv+0x2d7)[0xaa9467]
/home/projects/bbp00039/fastscape_build_files/aspect-22feb22-dealii-10.0-pre-fastscape-esther(_Z13run_simulatorILi2EEvRKNSt7__cxx1112basic_stringIcSt11char_traitsIcESaIcEEES7_bbb+0x279)[0x1572f89]
/home/projects/bbp00039/fastscape_build_files/aspect-22feb22-dealii-10.0-pre-fastscape-esther(main+0x715)[0x9f54b5]
/lib64/libc.so.6(__libc_start_main+0xf5)[0x2aaad161c555]
/home/projects/bbp00039/fastscape_build_files/aspect-22feb22-dealii-10.0-pre-fastscape-esther[0x9fc9b7] 

[Djneu]

Hi all, I've hopefully addressed all the previous comments and fixed the indent issue. I think the plugin is now in a place for a second round of reviews. I'll note, there is a Named_VTK.f90 file that I have not updated and is still a bit of a mess, I'm hoping this will eventually be included within the fastscape code so I've left it alone for now. @gassmoeller @tjhei @naliboff @anne-glerum

[gassmoeller]

Also you should definitely add a changelog entry for this, and give us a reference to any publications that should be cited if users enable the fastscape coupling (we will add them here: https://aspect.geodynamics.org/citing.html?ver=release&src=www)

[Djneu]

Hi all, thank you for all the suggestions! Other than a few open questions I think I've addressed all the previous comments and the code is ready for another review. @anne-glerum @naliboff @gassmoeller @tjhei @bangerth

[gassmoeller]

Oh, you also mentioned open questions in your last post. Can you summarize which questions are open? Most previous comments are marked as resolved or outdated, so it is a bit hard to keep track of what is resolved and what is still open.

[tjhei]

Instead of us bundling a patched fastscape or providing these complicated instructions, can you maybe fork fastscape into your own github repo and apply your patches?

[Djneu]

Hi all, thanks for all the suggestions! I've tried to address most of the new comments and updated the ghost nodes quite a bit. For the named_vtk file, I've removed this and done as @tjhei suggested and made the changes in the default branch here: https://github.com/Djneu/fastscapelib-fortran Also, note that at the moment "export GFORTRAN_CONVERT_UNIT='big_endian'" must be used for the VTK files to work properly.

There are a few open questions, the more major ones being:

1. Whether there is a better way to do the quadrature rule for having enough interpolation points to fill the fixed fastscape resolution mesh. https://github.com/geodynamics/aspect/pull/4396#discussion_r1121663580

2. Somewhat related to the last question, at the moment having precision be a user-defined variable and the variable "use_velocities" was a work around to be able to use fastscape with the free surface. This comes with its own issues with the mesh not being fixed using normal projection, and losing topography resolution when interpolating from a lower resolution ASPECT mesh onto a higher resolution FastScape one. Generally, we've always let fastscape handle advection, in which case should we keep this option available and possibly try to improve it at some point or remove it?

3. Right now the restart files are still written outside the general ASPECT checkpointing, and I was having trouble getting the pre_checkpoint function to work. https://github.com/geodynamics/aspect/pull/4396#discussion_r1124365806

Some minor questions that could be changed but are maybe not necessary: https://github.com/geodynamics/aspect/pull/4396#discussion_r1157104612 - mass flux. (Rene) https://github.com/geodynamics/aspect/pull/4396#discussion_r905067660 - boundary condition integers (John) https://github.com/geodynamics/aspect/pull/4396#discussion_r1124375375 - get time. (Anne) https://github.com/geodynamics/aspect/pull/4396#discussion_r1081231170 - years vs. seconds (Anne) https://github.com/geodynamics/aspect/pull/4396#discussion_r774690022 - assert throw for ghost nodes. (Anne) https://github.com/geodynamics/aspect/pull/4396#discussion_r774692317 - bottom row for table. (Anne) https://github.com/geodynamics/aspect/pull/4396#discussion_r774695522 - ghost node periodicity (Anne)

A final thing to note, the written FastScape VTK's do not include a time. There is a separate script we generally run after the model finishes to convert these from VTK to a VTS file that includes the correct timing to match ASPECT. I can include this later, but maybe this conversion is something we could also have ASPECT do automatically after FastScape writes the file?

[bangerth]

Separately, I think that this PR is currently based on some commit from 2022. Can you rebase to current main?

[anne-glerum]

@Djneu I pushed some documentation changes to your branch, this means you'll have to pull before you can push again.

[bangerth]

Thanks, @Djneu , for addressing my comments. Can you push your code changes? I've got a couple more questions where I wanted to look at how the new code looks that you mention in your replies.

When pushing, make sure you don't overwrite @anne-glerum 's commits.

[Djneu]

Hi @bangerth, thank you for all the comments! I've updated my branch to where I've gotten to so far, which should be most of the changes in fastscape.cc and fastscape.h that were mentioned.

I've also added a python script in the contrib folder, this is something that is necessary to visualize ASPECT and FastScape visualizations at the same time, as the general FastScape vtk files do not write a time. The script converts the VTK files into VTS using the solution.pvd times so they will match up.

[Djneu]

Hi all, Thanks again for all the suggestions! Other than a few minor comments (e.g., how to handle the GFORTRAN_CONVERT_UNIT), I think I've addressed all the comments and the code is ready for another review. @bangerth @gassmoeller @anne-glerum @naliboff @tjhei

[bangerth]

But it turns out that I can't merge, even though I want, github just doesn't let me :-) You'll have to address the following things:

• Run the ./contrib/utilities/indent script.
• There's a test that fails and that looks like this:

 302/1045 Test  #302: free_surface_VE_cylinder_2D_loading_fixed_elastic_dt_gmg ...........***Failed    0.66 sec
*** "Test free_surface_VE_cylinder_2D_loading_fixed_elastic_dt_gmg failed": ***
[1/1] Generating output-free_surface_VE_cylinder_2D_loading_fixed_elastic_dt_gmg/screen-output
FAILED: output-free_surface_VE_cylinder_2D_loading_fixed_elastic_dt_gmg/screen-output 
cd /__w/aspect/aspect/build/tests && /__w/aspect/aspect/tests/cmake/timeout.sh 600s mpirun -np 2 /__w/aspect/aspect/build/aspect /__w/aspect/aspect/build/tests/free_surface_VE_cylinder_2D_loading_fixed_elastic_dt_gmg.x.prm > /__w/aspect/aspect/build/tests/output-free_surface_VE_cylinder_2D_loading_fixed_elastic_dt_gmg/screen-output.tmp && /__w/aspect/aspect/tests/cmake/default screen-output </__w/aspect/aspect/build/tests/output-free_surface_VE_cylinder_2D_loading_fixed_elastic_dt_gmg/screen-output.tmp >/__w/aspect/aspect/build/tests/output-free_surface_VE_cylinder_2D_loading_fixed_elastic_dt_gmg/screen-output

--------------------------------------------------------
An error occurred in line <1565> of file </opt/tmp/unpack/deal.II-master/include/deal.II/lac/la_parallel_vector.h> in function
    Number dealii::LinearAlgebra::distributed::Vector<Number, MemorySpace>::operator()(dealii::LinearAlgebra::distributed::Vector<Number, MemorySpace>::size_type) const [with Number = double; MemorySpace = dealii::MemorySpace::Host; dealii::LinearAlgebra::distributed::Vector<Number, MemorySpace>::size_type = unsigned int]
The violated condition was: 
    partitioner->in_local_range(global_index) || vector_is_ghosted == true
Additional information: 
    You tried to read a ghost element of this vector, but it has not
    imported its ghost values.

Do you have an idea why that might fail? Several other tests fail with the same error. It is possible that if you rebase to the current main branch, the error will go away. Can you try a rebase?

• While you're there, and doing the rebase anyway, would you mind merging the commits into one (or one for the changes, one for tests, one for the changelog entry, one for documentation, etc -- i.e., into separate, self-contained commits)?

Thanks for your patience with this!

[Djneu]

Hi @bangerth, thank you for the reply! I've rebased the branch, merged all the commits into one, and ran the indent script. For the test error, after rebasing and running the tests I still have a failure on that test but it looks like its from differing outputs and not from the dealii error, so hopefully the rebase fixed it!

I've written up some documentation for using and installing the coupling (although I probably need to update it again). Should this be added in a separate PR?

[gassmoeller]

/rebuild

[Djneu]

Hi all, It looks like both the failed checks are related to the fastscape_eroding_box cookbook. It makes sense this won't run as fastscape.cc isn't compiled. Is there anything I should do for this?

[gassmoeller]

Ah, I forgot about this, yes, there are two simple fixes:

1. The cookbook markdown file was not listed in the cookbook index, and it didnt contain a proper markdown heading, this caused warnings in our documentation tester.
2. We automatically test if all cookbooks prms at least conform to proper .prm file syntax. Since fastscape is not compiled by default the cookbook you added had to manually be listed as exempt from that test.

I pushed a commit that should fix both issues. Let's see what the testers say, if they pass, I will merge.