[WIP] Add prescribed dilation to visco plastic #5803
Conversation
| set End time = 0e3 | ||
| set Use years in output instead of seconds = true | ||
| #set Nonlinear solver scheme = single Advection, iterated Newton Stokes | ||
| set Nonlinear solver scheme = single Advection, iterated Stokes |
There was a problem hiding this comment.
I tried Newton, and the convergence was not great, but will try again.
There was a problem hiding this comment.
Overall, the problem is very difficult to solve with prescribed dilation on. It would be good to compare convergence rates with @cedrict.
| @@ -295,8 +295,47 @@ namespace aspect | |||
| if (ElasticAdditionalOutputs<dim> *elastic_out = out.template get_additional_output<ElasticAdditionalOutputs<dim>>()) | |||
| { | |||
| elastic_out->elastic_shear_moduli[i] = average_elastic_shear_moduli[i]; | |||
| } | |||
| } | |||
There was a problem hiding this comment.
The most recent commits were indented.
naliboff
force-pushed
the
add_prescribed_dilation_to_visco_plastic
branch
from
495d7be
to
64d2166
Compare
naliboff
force-pushed
the
add_prescribed_dilation_to_visco_plastic
branch
from
64d2166
to
102c5ab
Compare
|
@tjhei - I've updated the PR to (1) remove the benchmark prm, (2) add a test case based off of the benchmark prm, (3) improve the code structure a bit where we calculate the dilation rate. The convergence behavior is still poor for the cases with prescribed dilation turned on, but my recollections is that @cedrict also ran into these issues as well with his models? For ASPECT, a minimum resolution (5 refinement levels for the test case) is required to get to 1e-3 for standard picard iterations (defect correction does equally poorly). The plot below shows the results for plastic dilation turned off (left column) versus turned on (right column) for different global refinement levels. The model with plastic dilation turned on and a refinement level of 4 has not yet converged. However, the plastic dilation is still clearly producing relatively resolution-independent shear band angles. At this stage, it is an experimental feature that will require additional (lots?) of work for robust usage. Do we want to merge this features in a piecemeal fashion, or all at once when we are confident it works robustly?
|
naliboff
force-pushed
the
add_prescribed_dilation_to_visco_plastic
branch
from
102c5ab
to
e7ec9c7
Compare
|
|
||
| /** | ||
| * The strain rate invariant | ||
| */ |
| @@ -297,6 +297,27 @@ namespace aspect | |||
| elastic_out->elastic_shear_moduli[i] = average_elastic_shear_moduli[i]; | |||
| } | |||
| } | |||
|
|
|||
| // Prescribed dilation | |||
| double average_friction_angle = 0.; | |||
|
|
||
| if (prescribed_dilation) | ||
| { | ||
| average_friction_angle = MaterialUtilities::average_value(volume_fractions, |
There was a problem hiding this comment.
I think you only need to compute this if you are yielding
| MaterialModel::PrescribedPlasticDilation<dim> | ||
| *prescribed_dilation = out.template get_additional_output<MaterialModel::PrescribedPlasticDilation<dim>>(); | ||
|
|
||
| if (prescribed_dilation) |
There was a problem hiding this comment.
should this be an error if you enable dilation but you don't have the outputs? Or does this break somewhere?
There was a problem hiding this comment.
My thought was we would simply not fill the outputs if they don't exist with this logic. Is there a better way to ensure this is the case?
| if (plastic_yielding) | ||
| dilation_rate = 2. * std::sin(average_friction_angle) * isostrain_viscosities.strain_rate_invariant; | ||
|
|
||
| prescribed_dilation->dilation[i] = dilation_rate; |
There was a problem hiding this comment.
Why don't you assign it directly here?
There was a problem hiding this comment.
Also, do you need to set it =0; otherwise?
There was a problem hiding this comment.
Why don't you assign it directly here?
The default value is set to 0 on line 303, and it is not assigned directly here as I was not sure how the logic for plastic_yielding = true or false would be worked in on a single line. Open to suggestions.
Supersedes #3119, but is only open now for discussion and testing.
@cedrict @tjhei @bobmyhill - The PR contains a working example (it runs).
The code can be cleaned up quite a bit (see commented ToDo), and one question to consider is whether we want to calculation the dilation rate from an average friction angle across compositions (current approach) or calculate a dilation rate for each composition and then average them (would require refactoring).
For some reason the dilation values are not being written out as a visualization field (they should be when selecting
named additional outputs), but I probably just forgot to enable something.Time to start testing :)
For all pull requests:
For new features/models or changes of existing features: