Need the ability for UserObjects to contribute to the Jacobian

[dschwen]

• Add flag to Assembly that is set to true as soon as even one UO requires the shapes (we have a setter that can only set to true - there is no reason to ever unset this).
• Query that flag in ComputeUserObjectsThread (we'll retrieve a const ref from assembly)
• call FEProblem::prepareShapes in ComputeUserObjectsThread::onElement if the flag is true
• The new base class ShapeElementUserObject will provide references to _phi and _grad_phi

Ping @SudiptaBiswas

[dschwen]

Slight change of plans. The ShapeElementUserObjects will register each of their coupled variables using a function in assembly. Assembly will maintain a vector of variable IDs for which shape functions are requested. ComputeUserObjectsThread will retrieve a const ref to that vector and loop over it in ::onElement calling FEProblem::prepareShapes for each entry.

It is not pretty, so I'd welcome better ideas!

[permcody]

This does sound expensive. As long as the normal code path isn't impacted,
this should be OK.

On Tue, Nov 3, 2015 at 10:24 AM Daniel Schwen notifications@github.com
wrote:

Slight change of plans. The ShapeElementUserObjects will register each of
their coupled variables using a function in assembly. Assembly will
maintain a vector of variable IDs for which shape functions are requested.
ComputeUserObjectsThread will retrieve a const ref to that vector and
loop over it in ::onElement calling FEProblem::prepareShapes for each
entry.

It is not pretty, so I'd welcome better ideas!

—
Reply to this email directly or view it on GitHub
#5913 (comment).

[dschwen]

I'll push a commit which implements this, now. This needs to be in Assembly. Otherwise the UO constructors from all threads will cause race conditions.. unless I register only for _tid == 0...

[dschwen]

I'm protecting the loop over the vector entries by a check against _fe_problem.currentlyComputingJacobian(). That way we only add that overhead once per non-linear iteration.

[friedmud]

Guys: please hold off on this! We need to think about this and design it properly.

This is not a thing to just "DO". It's a thing that must be done carefully. Please create a Ticket for what you NEED (not how to do it) and let's discuss possibilities. An error in design here will cost us a lot of time in maintenance down the road...

[permcody]

We won't merge this until we have a chance to talk with you. Although, we probably should have that conversation fairly soon.

[friedmud]

It's not about merging. Code shouldn't be written until we've had a design
discussion. Once code is written it's hard to change course.

UserObjects are not designed to contribute to Jacobians. If they are going
to start doing that, that is going to be a MAJOR change.

Firstly, we need to think about whether or not we want to allow it all. And
if we do then we need to talk about what the real requirements are, what
deliverable they're tied to and when this capability is due.

Next, we need to double check that we should be doing this in a
UserObject... or whether we should create a new System.

Third, we need to make sure that there is design coherency between how
Jacobians are computed in Kernels, BCs, DiracKernels and UserObjects. There
should be just one design pattern that users need to learn.

Finally, we need to think about alternative uses, how it's going to work
with nodal, and general user objects, etc.

There is a LOT to talk about before code gets written.

Derek

On Tue, Nov 3, 2015 at 4:36 PM Cody Permann notifications@github.com
wrote:

We won't merge this until we have a chance to talk with you. Although, we
probably should have that conversation fairly soon.

—
Reply to this email directly or view it on GitHub
#5913 (comment).

[permcody]

Alright, fair enough, @dschwen is already head deep into this.

Let's talk about alternatives, he needs to contribute to the Jacobian, but in a global sense. Perhaps there's an easy way to adapt our existing systems such that he could accumulate a value in a Kernel and then add it in at the end on a "post-Jacobian" call-back or something. I'm sure Daniel is open for ideas.

[dschwen]

Yeah, I'm in way deep :-). So, short story is, I didn't think this would be a high impact change. We're adding a little bit of flexibility enabling a specific UserObject type that can calculate derivatives of some quantities it computes with respect to the DOFs. This data can be used by a kernel that uses the UO quantity to construct the correct Jacobian entries.

So what do we need (as opposed to how). We have a PDE (an advection term) for which one parameter (advection velocity) depends on an expression that contains an integral over the simulation domain. The integral is over a force density, which peaks wherever certain particles in the phase field simulation touch. I.e. the contributions to the integral are very heterogenous over the simulation domain. We believe that this is the reason we really need the Jacobian contributions from that integral (replacing the integral with a constant advection velocity makes the problem converge as opposed to being unsolvable for example).

Integral made us think Postprocessor (we actually need a UO as we are integrating vector quantities). The ElementUserObject could assemble a vector of all Jacobian contributions that the Advection Kernel could possible need with very little overhead (no additional looping over elements!).

I'd be really happy if we can come up with a more elegant design on how to implement this. But keep in mind, that the UserObject won't be writing to the Jacobian matrix, it will need to provide terms that a Kernel then uses. UserObject interfaces are already completely user defined, which means the Jacobian term interface must be completely user defined as well.

[YaqiWang]

@dschwen As a reminder, have you tried FD solve_type to see if the exact Jacobian assembled through PETSc finite difference helps? We had similar thoughts (postprocessors are coupled into kernels to take over constant values) before but it turns out the other part of Jacobian may not be that important.

[SudiptaBiswas]

I have tried FD before, while developing this model. That did not help in solve/convergence. Will recheck how this recent development impacts.

[dschwen]

Wait... if FD does not help then fully correct jacobians will not help either!

[SudiptaBiswas]

Let me re-check it.

[SudiptaBiswas]

Sorry, I was wrong. I did not try solve-type FD, I have tried FDP preconditioner with solve-type = Newton & pc-type lu. Will see how FD works.

[SudiptaBiswas]

Okay, here it is (with existing RigidBodyMotion UO & kernels).

1. solve-type = FD with FDP preconditioner converges fine.
2. solve-type = PJFNK, FDP preconditioner, converges fine with small time step
3. solve-type = NEWTON, FDP preconditioner, takes more no. of non-linear iterations to converge.
4. Jacobian test results;
   Kernel for variable 'c':
   (0,1) Off-diagonal Jacobian for variable 'w' is slightly off (by 0.106266 %)
   Kernel for variable 'w':
   (1,0) Off-diagonal Jacobian for variable 'c' is slightly off (by 0.683286 %)
   (1,1) On-diagonal Jacobian is wrong (off by 62.6 %)

That being said, I think we are heading towards right direction.

[YaqiWang]

That is good! With solve_type = FD, you will not need FDP. Solve_type= PJFNK and full=true in FDP, should be equal to solve_type = FD. Does NEWTON+FDP converge to the same solution as PJFNK+FDP? If FDP full=true, they should.

[SudiptaBiswas]

I was expecting the same, however solve_type = FD with no preconditioner specified was throwing some Petsc error. I did not dig into this much as it went off with addition of preconditioner.
Solve_type= PJFNK and full=true in FDP is what the 2nd case of my previous comment is. Newton converges to similar order of residual norm after cutting down the time step & uses more no. of nonlinear iterations.

[lindsayad]

Shouldn't this issue be closed per #7216?

[SudiptaBiswas]

Yes, this should be closed.

[permcody]

Closed by #7216