On Mon, Aug 28, 2017 at 5:15 AM, Maximilian Hartig <imilian.har...@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.har...@gmail.com> > wrote: > > > On 17. Aug 2017, at 13:51, Matthew Knepley <knep...@gmail.com> wrote: > > On Thu, Aug 17, 2017 at 3:13 AM, Maximilian Hartig <imilian.hartig@gmail. > com> wrote: > >> Thanks for your help Matt. >> >> On 16. Aug 2017, at 16:17, Matthew Knepley <knep...@gmail.com> wrote: >> >> On Wed, Aug 16, 2017 at 8:56 AM, Maximilian Hartig < >> imilian.har...@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/