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

Maximilian Hartig imilian.hartig at gmail.com
Mon Aug 28 09:51:04 CDT 2017 

> On 28. Aug 2017, at 12:31, Matthew Knepley <knepley at gmail.com> wrote:
> 
> On Mon, Aug 28, 2017 at 5:15 AM, Maximilian Hartig <imilian.hartig at gmail.com <mailto: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 <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.

ok, so now I have a tspoststep function in which I do:

TSGetDM(ts, &dm);
PetscObjcetQuery((PetscObject) dm, “A”, (PetscObject*) &copied_solution);
TSGetSolution(ts, &solution);
VecCopy(solution, copied_solution);

but Petsc complains that the vectors have different sizes. I guess this is due to the Dirichlet BC I imposed. I tried imposing the same boundary condition on the auxiliary field which made the vectors the same length. However it seemed to mess up the way auxiliary fields get accessed. Any idea on how I can get the “full” solution vector where the boundary values are filled in? I tried writing it to a hdf5- file and then calling VecRead on it but I cannot get the hdf5 file format to work.
>  
> 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 <http://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/TS/TSSetPreStep.html#TSSetPreStep>
That seems very useful as well. Are there any examples for the sub solve?

Thanks,
Max
> 
> 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 <mailto:imilian.hartig at gmail.com>> wrote:
>> 
>>> 
>>> On 17. Aug 2017, at 13:51, Matthew Knepley <knepley at gmail.com <mailto: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/>
> 
> 
> 
> -- 
> 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/20170828/382483e2/attachment-0001.html>

More information about the petsc-users mailing list