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

[Maximilian Hartig]

> 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 <mailto:imilian.hartig at gmail.com>> wrote:
> Thanks for your help Matt.
>> On 16. Aug 2017, at 16:17, Matthew Knepley <knepley at gmail.com <mailto:knepley at gmail.com>> wrote:
>> 
>> On Wed, Aug 16, 2017 at 8:56 AM, Maximilian Hartig <imilian.hartig at gmail.com <mailto: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/ <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/ <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/ <http://www.caam.rice.edu/~mk51/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20170822/4a8f3b55/attachment-0001.html>