On Mon, Dec 2, 2013 at 4:29 PM, Jed Brown <[email protected]> wrote:
> Geoffrey Irving <[email protected]> writes: > > > For weak form PDEs which arise as the gradient of an integrated energy > > function, it would be great to have a routine to perform the > > appropriate integral. It would be best to structure this as a generic > > "perform an integral" function, evaluating > > > > E(f) = sum_i w_i f(u(x_i), grad u(x_i)) > > There is a SNES interface for quadratic functionals, SNESSetObjective(), > and it has been plumbed down to DMDASNESSetObjectiveLocal(), but there > is no DMSNESSetObjectiveLocal. I would add that, and a > DMPlexComputeResidualFEM(). > > Actually, I would have made DMPlexSNESSetFunctionFEM() and now > DMPlexSNESSetObjectiveFEM(), which configures the SNES appropriately, > rather than calling SNES methods and passing in pointers to library > functions. All the other DMs work this way. Matt, why did you do it > differently with DMPlex? > When you and Peter changed the way SNES used to work, I asked you how I should change Plex. You told me to do this: ierr = DMSNESSetFunctionLocal(dm, (PetscErrorCode (*)(DM,Vec,Vec,void*)) DMPlexComputeResidualFEM, &user);CHKERRQ(ierr); ierr = DMSNESSetJacobianLocal(dm, (PetscErrorCode (*)(DM,Vec,Mat,Mat,MatStructure*,void*)) DMPlexComputeJacobianFEM, &user);CHKERRQ(ierr); which I did. Moreover, this is how it always worked, namely specifying the local function. What exactly are you proposing to replace this? > > > for a set of FE fields u (plus auxiliary fields) and some quadrature > > rule. The quadrature rule might need to be passed in independently of > > the PetscFE objects, since the integral might tie together several > > different fields connected to several different PetscFE's, conceivably > > with distinct quadrature rules. The integral is only valid as a > > discrete energy if all the quadrature rules match, but there might be > > other integral application where this doesn't hold. > > How is any of this different from evaluating residuals? Each element > computes a scalar now, all of which are summed together, but otherwise I > think it's the same mechanics. It looks the same as DMPlexComputeL2Diff(). > > Relatedly, how does one ensure that a PetscFE object is Galerkin > > (matching primal and dual spaces) so that the energy to be discretely > > valid? > > How does one ensure that the weak form itself is symmetric and how do > you ensure that it is not indefinite (like a saddle; the Lagrangian is > not the energy)? We can do consistency checks, like verifying that the > provided residual function is the gradient of the objective (energy). > We do not have any way to check that the spaces are the same, and I think its not trivial. Here I think we fall back on the user giving us the correct input. Matt -- 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
