Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Free surface tractions #3058

Merged
merged 2 commits into from Feb 6, 2020

Conversation

fionaclerc
Copy link
Contributor

This pull request add an exception to allow boundary tractions to be applied on the free surface. It also includes benchmarks for the surface deformation of viscous and viscoelastic materials in response to surface loading.

Copy link
Contributor

@naliboff naliboff left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Outstanding, thank you for contributing these benchmarks! I've done an initial review and on the whole this is PR is in good shape.

A few general comments and requests.

  1. Please remove the example_1.prm
  2. Would you mind include a python script in the benchmarks folder that calculates the analytical solution?

Nice work and great to see the viscoelasticity implementation is working at intended!

solution found in the 'NL82_soln' directory, using the provided
gnuplot script ('compare_VE_def' for maximum surface deflection
through time, 'compare_VE_def_profile' for deflection of profile
through time).
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would you be willing to put a python script in this directory that calculates the analytical solution and writes out these files?

@@ -0,0 +1,7 @@
set title "Depression of surface in response to loading (inst. from 0 to 1000 years), over viscoelastic medium"
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you add a brief comment here explaining what this file does?

@@ -0,0 +1,11 @@
set title "Depression of surface in response to loading (inst. from 0 to 1000 years), over viscoelastic medium"
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same as above.


# Mesh refinement specifications
subsection Mesh refinement
set Initial adaptive refinement = 1
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do the results change significantly if you do not do the adaptive refinement? If so, it might be worth putting in a brief comment mentioning this.

# r0 is load radius, H0 is load height, t1 is time load is removed,
# rhoi is density of ice/load, g is gravitational acceleration
set Function expression = 0; if (x<r0 && t<t1,-g*H0*rhoi*(1-t/t1), 0)
#set Function expression = 0; if (x<r0 ,if(t<t0,-g*H0*rhoi,if(t<t1,-g*H0*rhoi*(1-(t-t0)/t1),0)), 0)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove this line.

# Temperature boundary conditions
subsection Boundary temperature model
set Fixed temperature boundary indicators = bottom, top, left, right
set List of model names = box
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You can change this line to set List of model names = initial temperature, which would allow deleting the Box subsection below.

# Material model
subsection Material model

set Model name = simple
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For simplicity, you could use the viscoelastic material model here and set the Elastic shear moduli to something like 1e50, which for all intensive purposes produces viscous behavior (e.g., negligible viscoelastic contribution).

end
end

# Number and name of compositional fields
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Unless you use the Viscoelastic material model below, there is no need to specify these compositional fields.

@@ -2040,7 +2040,7 @@ namespace aspect
std::set<types::boundary_id> boundary_indicator_lists[6]
= { boundary_velocity_manager.get_zero_boundary_velocity_indicators(),
boundary_velocity_manager.get_tangential_boundary_velocity_indicators(),
parameters.free_surface_boundary_indicators,
// parameters.free_surface_boundary_indicators,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please delete this line.

@naliboff
Copy link
Contributor

/rebuild

@gassmoeller
Copy link
Member

@fionaclerc do you want to talk about how to address the comments?

@naliboff
Copy link
Contributor

@fionaclerc - is this ready for another review?

Copy link
Contributor

@naliboff naliboff left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for updating and reorganizing all of this @fionaclerc! I think this is very close to being ready.

After you address the comments here, I would like to take one last look over to check for minor spelling errors, etc. Overall, the PR is in great shape and will be a really useful contribution to the code.

After merging, I think a cookbook on viscoelasticity is in order.

General comments and questions:

  1. Does example_1.prm in the main folder still need to be included? If so, can you add some brief documentation (and perhaps rename)? If not, go ahead and delete this file.

  2. Please add a changelog entry (doc/modules/changes) documenting the changes to the source code and benchmarks!

Changing the stress instantaneously as the load is applied/removed
leads to oscillations in the elastic response. This can be addressed
by setting a fixed elastic timestep larger than the numerical
timestep and turning on stress averaging:
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do the oscillations still appear if you use a small CFL value and a Maximum time step that is equal to the elastic time step? In other problems, that scenario has produced reasonable results.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Using a smaller CFL value does not reduce the oscillations. Setting the maximum time step equal to the elastic time step does work, as in the parameter files “free_surface_VE_cylinder_*D_loading_fixed_elastic_dt.prm”. I am correcting/clarifying the README file.

r0 is the load radius), for t>0. The load is fully removed by t=t1,
and the surface relaxes. This is done both in a 2-D and 3-D geometry
(by symmetry the load is centered on the left boundary or left/front
corner). The input files are:
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you provide a bit more detail in this paragraph about the different loading/unloading scenarios (e.g., instantaneous verse linear)?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes.

@@ -0,0 +1,35 @@
Benchmark of infinite viscous half-space loaded/unloaded by an
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Consider adding a note here and in the parameter files that the viscous benchmark is identical to the viscoelastic benchmarks, with the exception of the rheology (e.g., viscous verse viscoelastic).

@fionaclerc
Copy link
Contributor Author

I have updated the README files and some comments to address the review. I also added a note in the README files about the mismatch between the analytical and numerical solutions at the far (right) boundary due to the free-slip condition.

@fionaclerc fionaclerc closed this Dec 19, 2019
@fionaclerc fionaclerc reopened this Dec 19, 2019
@tjhei
Copy link
Member

tjhei commented Dec 20, 2019

thanks for continuing the work on this! Sadly, there have been edits in free_surface.cc merged since you started your branch. You will need to rebase your changes to the current master. Note that the free_surface.cc got moved into a different folder, so you likely have to fix the conflicts manually. Let us know if you need help.

@naliboff
Copy link
Contributor

@fionaclerc and I are working on fixing the rebase, will rebuild and test after that is done.

@fionaclerc fionaclerc force-pushed the free_surface_tractions branch 2 times, most recently from 9dbcf15 to 9dc0599 Compare January 30, 2020 22:03
@naliboff
Copy link
Contributor

naliboff commented Feb 4, 2020

@fionaclerc - This looks great, well done and thanks for sticking through all of the revisions! These benchmarks will be incredibly useful for a large group of users.

The last step is to do a final commit adding a changelog entry in the doc/modules/changes directory. A good example for this type of addition is: doc/modules/changes/20190529_euen.

Let me know if you have any questions about this and as soon the change log file is added, this will be ready to merge.

Co-Authored-By: naliboff <john.naliboff@gmail.com>
@naliboff naliboff merged commit ef542ec into geodynamics:master Feb 6, 2020
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

4 participants