[petsc-users] accessing fields from previous, converged solution.
Matthew Knepley
knepley at gmail.com
Mon Aug 28 16:04:04 CDT 2017
On Mon, Aug 28, 2017 at 10:51 AM, Maximilian Hartig <
imilian.hartig at gmail.com> wrote:
>
> 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> 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.
>
>
> 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.
>
Okay, this should be better documented. The auxiliary fields are all local
vectors. This makes it easy to project
them to any cell. It means that above you will replace VecCopy with
DMGlobalToLocalBegin(dm, solution, copied_solution);
DMGlobalToLocalEnd(dm, solution, copied_solution);
and you will need a local vector for copied_solution. Sorry bout that.
> 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
>
> That seems very useful as well. Are there any examples for the sub solve?
>
There is one in PyLith (dynamic fault friction), but I cannot think of one
in the PETSc examples right now.
Thanks,
Matt
> 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>
>> 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/
>
>
>
--
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/d53f1b00/attachment.html>
More information about the petsc-users
mailing list