[petsc-users] accessing fields from previous, converged solution.

[Matthew Knepley]

On Mon, Aug 28, 2017 at 5:15 AM, Maximilian Hartig <imilian.hartig at gmail.com
> wrote:

> Sorry to bother you again, but I have checked and rechecked my formulation
> and it corresponds to what I have found in literature. The only difference
> is that my flow criteria is not evaluated using the trial stress and the
> current yield stress but rather the corrected stress and the updated yield
> stress. I am convinced this is what causes the problem.
> Unfortunately I do not have any idea on how to go about implementing this.
> First I don’t know how to access the solution from the past iteration
>

I think you can do what you want using


http://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/TS/TSSetPostStep.html

In your function, you would TSGetSolution() and then copy it into an
auxiliary vector you
set. This avoids the problem of composition you had before.


> and secondly I don’t know how to do the elastic predictor step.
>

You can always create another SNES and do a subsolve, maybe in


http://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/TS/TSSetPreStep.html#TSSetPreStep

which you then copy into an auxiliary vector, just as above.

Tell me if I am missing something.

  Thanks,

     Matt


> I tried doing it with auxiliary fields that get updated every timestep but
> I gather that is not a good idea. Could you please point me in some
> direction?
>
> I appreciate you patience with my questioning.
>
> Thanks,
> Max
>
> On 22. Aug 2017, at 13:52, Maximilian Hartig <imilian.hartig at gmail.com>
> wrote:
>
>
> On 17. Aug 2017, at 13:51, Matthew Knepley <knepley at gmail.com> wrote:
>
> On Thu, Aug 17, 2017 at 3:13 AM, Maximilian Hartig <imilian.hartig at gmail.
> com> wrote:
>
>> Thanks for your help Matt.
>>
>> On 16. Aug 2017, at 16:17, Matthew Knepley <knepley at gmail.com> wrote:
>>
>> On Wed, Aug 16, 2017 at 8:56 AM, Maximilian Hartig <
>> imilian.hartig at gmail.com> wrote:
>>
>>> Hello,
>>>
>>> I have a problem with several fields that I solve with PetscFE and TS. I
>>> now need to access the solution from the previous timestep to compute the
>>> residual for the current timestep.
>>> I tried a TSMonitor with the following code in it:
>>>
>>> TSGetDM(ts,&dm);
>>> DMClone(dm,&dm_aux);
>>> DMGetDS(dm,&prob_aux);
>>> DMSetDS(dm_aux,prob_aux);
>>> DMCreateGlobalVector(dm_aux,&old_solution);
>>> VecCopy(u,oldsolution);
>>> PetscObjectCompose((PetscObject) dm, “A”, (PetscObject) old_solution);
>>>
>>> VecDestroy(&old_solution);
>>> DMDestroy(&dm_aux);
>>>
>>> hoping that it would create an auxiliary field that I could access in
>>> the evaluation of the residual. It did that but messed with the
>>> discretisation of the initial problem in some way. So I figure that adding
>>> auxiliary fields to a dm after having fed it to a TS context is not
>>> something you should be doing.
>>>
>>> Is there a way to access the fields of the solution for the previous
>>> timestep during the evaluation of the current residual?
>>>
>>
>> First, I can show you how to do what you are asking for. I think you can
>> simply
>>
>> PetscObjectQuery((PetscObject) dm, "old_solution", (PetscObject *)
>> &old_solution);
>> if (!old_solution) {
>>   DMCreateGlobalVector(dm, &old_solution);
>>   PetscObjectCompose((PetscObject) dm, "old_solution", old_solution);
>> }
>> VecCopy(u, oldsolution);
>>
>>
>> Unfortunately, this produces an error and tells me “An object cannot be
>> composed with an object that was composed with it."
>>
>>
>> Second, I think a better way to do this than composition is to use
>>
>>   DMGetNamedGlobalVector(dm, "old_solution", &old_solution);
>>
>>
>> Yes, this seems to work and I do not get an error while running the
>> monitor. Also the discretisation seems to be fine. But the old solution now
>> is not available as an auxiliary field.
>>
>>
>> Third, I can say that I am profoundly troubled by this. This interferes
>> with the operation of
>> the time integrator, and I cannot see a reason for this. If you are
>> keeping track of the integral
>> of some quantity, I would update that in TSPostStep() and request that
>> integral instead of the
>> previous solution. We do this for some non-Newtonian rheologies.
>>
>>
>> I have an “if” condition in my residual which I need to evaluate to
>> determine the correct formulation for my residuals. I use actual fields to
>> keep track of the integrals. But when I use the actual field to evaluate
>> the “if” condition, due to the implicit nature of the solver I jump “back
>> and forth” over that condition. I need the “if” condition to be independent
>> from the field for which it determines the residual.
>> The actual application is as follows:
>>  I have, amongst others, a displacement and a stress field.
>> I evaluate for my stress given a displacement increment delta u.
>> If the resulting stress is > yield stress, I need a plasticity
>> formulation, if it is smaller, I can use elasticity.
>>
>
> Hmm, I think this is a formulation problem. You should be using the
> "correct" formulation at the implicit step, since otherwise
> the residual will be inconsistent with what you tested when evaluated at
> the new step. I have the same kind of elasto-plastic
> system in PyLith (http://www.geodynamics.org/cig/software/pylith/) which
> we solve implicitly without much problem.
>
>
> Possible, but I rechecked the formulation and it looks fine to me. In
> literature, this kind of algorithm uses an elastic predictor step to
> determine whether to use the elastic or the plastic formulation. The
> elastic “trial” stress is compared with the yield stress from the past
> timestep. In my current setup, I do not compare the elastic “trial” stress
> with the past yield stress but the current stress (including plastic
> correction) with the current yield stress. I think this is what leads to my
> problem. When I cheat and just tell petsc to use the plastic formulation
> from the time on where I expect plastic behaviour, I get convergence and
> the results look reasonable.
> So I tried to implement a flow condition which would not change during the
> nonlinear solver iteration steps. If the first, elastic predictor step
> results in yielding then, in my understanding, from there on we should
> expect plastic behaviour and not jump “back and forth” between plastic and
> elastic behaviour.
> I have taken a look at pylith and it seems to use an explicit formulation
> for dynamic plasticity. I try to use an implicit formulation solving for
> displacement, stress, plastic strain and the plastic consistency operator
> at once. I am not aware of a possibility within the petsc framework to have
> an internal solver to compute stress state and internal variables depending
> on the displacement rate and then return to the outer iteration to solve
> for the correct displacements. This is how I found it done in literature.
> Is there a possibility to have such an “internal” snes?
>
> Thanks,
> Max
>
>
>   Thanks,
>
>     Matt
>
>
>> Thanks,
>> Max
>>
>>
>>   Thanks,
>>
>>      Matt
>>
>>
>>> Thanks,
>>> Max
>>
>>
>>
>>
>> --
>> What most experimenters take for granted before they begin their
>> experiments is infinitely more interesting than any results to which their
>> experiments lead.
>> -- Norbert Wiener
>>
>> http://www.caam.rice.edu/~mk51/
>>
>>
>>
>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
> http://www.caam.rice.edu/~mk51/
>
>
>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

http://www.caam.rice.edu/~mk51/
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20170828/2ddb3314/attachment.html>