On Thu, Aug 17, 2017 at 3:13 AM, Maximilian Hartig <imilian.har...@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.hartig@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. 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/