_du_dot_du[_qp] in TimeDerivative::computeQpJacobian() might not be so long as the mass matrix supposed.

[lw4992]

Real TimeDerivative::computeQpJacobian() { return _test[_i][_qp]phi[_j][_qp]du_dot_du[_qp]; } It was designed to compute mass matrix, _du_dot_du[_qp] shoud be dt for BDF1 simply. My test is modified as moose/test/tests/dgkernels/2d_diffusion_dg/2d_diffusion_dg.i. I add a TimeDerivative kernel to test transient Jacobian using TimeDerivative::computeQpJacobian() like this [./time] type = TimeDerivative variable = u [../] But, when I print _du_dot_du[_qp] out, they were NOT my input dt. The reason maybe lie in moose/framework/src/base/MooseVariable.C. In line 563, 780 and 922, _du_dot_du[qp] += phi_local * du_dot_du(idx) should be _du_dot_du[qp] = du_dot_du(idx);

Please investigate this issue, and check other similar issue.

[friedmud]

It's been a while since anyone has thought about this stuff.

@andrsd can you take a look at this?

[andrsd]

Re = \partial u \over \partial t * phi_i = u^(n+1) - u^(n) / dt * phi_i. Let du = u^{n+1} - u^{n} / dt, then du_dot_du = phi_j/dt.

So, you are right that are multiplying by an extra phi_j, but du_dot_du should be 1/ dt.

[andrsd]

I made a patch that removes the extra phi_j multiplication and 25 tests in moose failed. I did not try INL's internal apps yet. I'd like to postpone this until next week when the HPC infrastructure is back online...

[YaqiWang]

One suggestion: maybe you want to print out the Jacobian using moose jacobian evaluation and petsc evaluation to make sure the phi_j should be removed.

[friedmud]

I agree with @YaqiWang - let's make sure this is the right thing to do by verifying with PETSc. I cannot imagine we have had this wrong for so long - it would have shown up all over the place....

[lw4992]

@YaqiWang- In fact, I found this issue when comparing with the Jacobian using moose jacobian evaluation and petsc evaluation. For steady case, they are same, but transient case DIFFERENT! It makes me a feel that multiplying an extra phi_j in TimeDerivative.C and MooseVariable.C. Then I modified the code as previous, it works correctly.

[andrsd]

So: (1) the patch that I did yesterday was wrong, so disregard the message about 25 test failing test. (2) I build a simple problem: du\over dt = f, f = 1 and u = t on d\Omega (exact solution is u = t, which BDF1 solves exactly with first order lagrange elements). I ran that with -pc_type lu (2a) Your "fix" computes that exactly:

 0 Nonlinear |R| = 6.388271e-01
    0 KSP Residual norm 6.388270501474e-01 
    1 KSP Residual norm 1.132958715776e-15 
 1 Nonlinear |R| = 1.064859e-15

1.77973e-10 = ||J - Jfd||//J|| 1.12845e-09  = ||J - Jfd||

(2b) The current MOOSE code computes that exactly:

 0 Nonlinear |R| = 6.388271e-01
    0 KSP Residual norm 6.388270501474e-01 
    1 KSP Residual norm 1.132958715776e-15 
 1 Nonlinear |R| = 1.064859e-15

1.77973e-10 = ||J - Jfd||//J|| 1.12845e-09  = ||J - Jfd||

This shows the convergence and the last line is the verification with PETSc finite differencing, solve_type = 'NEWTON' in both cases.

I think that you got confused by not noticing the double loop that interpolates the du_dot_du vector into quadrature points. At least, I did not notice yesterday and that's why I thought we were multiplying by an extra phi. Well, we are not.

So, I'd like to see the input file you used that lead to the conclusion our du_dot_du is wrong...

[lw4992]

Only residual printing out can not indicate the issue since the Jacobian matrix size is so small. The residual convergence are almost same even though matrix is different. You should display the Jacbian matrix both hand-coded and FD with petsc by

 petsc_options = '-ksp_monitor -ksp_view -snes_test_display'
 petsc_options_iname = '-pc_type -snes_type'
 petsc_options_value = 'lu test'  

Here are my tests. The first case follows moose/test/tests/dgkernels/2d_diffusion_dg/2d_diffusion_dg.i. It is a steay case. Add the following block to display Jacobian matrix.

[Preconditioning]
    [./SMP]
      type = SMP
      full = true
          petsc_options = '-ksp_monitor -ksp_view -snes_test_display'
          petsc_options_iname = '-pc_type -snes_type'
          petsc_options_value = 'lu test'
    [../]
[]

It's OK!

Norm of matrix ratio 4.64326e-10 difference 3.2405e-08 (user-defined state)

Secondly, A TimeDerivative kernel was added into the input file to test TimeDerivative::computeQpJacobian() for transient case. The difference of matrix can be observed clearly.

Norm of matrix ratio 0.00477622 difference 0.333333 (user-defined state)

Lastly, with my "fix", test the transient case again. It's OK!

I don't know how to upload a file. Here is the copy of transient test.

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 2
  ny = 2
#  xmin = -1
#  xmax = 1
#  ymin = -1
#  ymax = 1
  xmin = 0
  xmax = 1
  ymin = 0
  ymax = 1
  elem_type = QUAD4
[]

[Variables]
  active = 'u'

  [./u]
    order = FIRST
    family = MONOMIAL

    [./InitialCondition]
      type = ConstantIC
      value = 1
    [../]
  [../]
[]

[Functions]
  active = 'forcing_fn exact_fn'

  [./forcing_fn]
    type = ParsedFunction
#    function = -4.0+(x*x)+(y*y)
#    function = x
#    function = (x*x)-2.0
    value = 2*pow(e,-x-(y*y))*(1-2*y*y)
#    function = (x*x*x)-6.0*x
  [../]

  [./exact_fn]
    type = ParsedGradFunction
#    function = x
#    grad_x = 1
#    grad_y = 0

#    function = (x*x)+(y*y)
#    grad_x = 2*x
#    grad_y = 2*y

#    function = (x*x)
#    grad_x = 2*x
#    grad_y = 0

    value = pow(e,-x-(y*y))
    grad_x = -pow(e,-x-(y*y))
    grad_y = -2*y*pow(e,-x-(y*y))

#    function = (x*x*x)
#    grad_x = 3*x*x
#    grad_y = 0
  [../]
[]

[Kernels]
  active = 'time diff abs forcing'

  [./time]
    type = TimeDerivative
    variable = u
  [../]

  [./diff]
    type = Diffusion
    variable = u
  [../]

  [./abs]          # u * v
    type = Reaction
    variable = u
  [../]

  [./forcing]
    type = UserForcingFunction
    variable = u
    function = forcing_fn
  [../]
[]

[DGKernels]
  active = 'dg_diff'

  [./dg_diff]
    type = DGDiffusion
    variable = u
    epsilon = -1
    sigma = 6
  [../]
[]

[BCs]
  active = 'all'

  [./all]
    type = DGFunctionDiffusionDirichletBC
    variable = u
    boundary = '0 1 2 3'
    function = exact_fn
    epsilon = -1
    sigma = 6
  [../]
[]

[Preconditioning]
    [./SMP]
      type = SMP
          full = true

    petsc_options = '-ksp_monitor -ksp_view -snes_test_display'
        petsc_options_iname = '-pc_type -snes_type'
    petsc_options_value = 'lu test'
        #petsc_options_iname = '-pc_type -pc_hypre_type'
    #petsc_options_value = 'hypre boomeramg'
    [../]

[]
[Executioner]
  type = Transient
  # Preconditioned JFNK (default)
  solve_type = 'NEWTON'

  [./Adaptivity]
    steps = 2
    refine_fraction = 1.0
    coarsen_fraction = 0
    max_h_level = 8
  [../]

  nl_rel_tol = 1e-10

[]

[Postprocessors]
  active = 'h dofs l2_err'

  [./h]
    type = AverageElementSize
    variable = u
  [../]

  [./dofs]
    type = NumDOFs
  [../]

  [./l2_err]
    type = ElementL2Error
    variable = u
    function = exact_fn
  [../]
[]

[Outputs]
  file_base = out
  output_initial = true
  exodus = true
  csv = true
[]

[lw4992]

Moreover, if print out _du_dot_du at quadrature points, some of them were not 1/dt (default = 1)

Real TimeDerivative::computeQpJacobian()
{
  std::cout << "_qp: "<<_qp <<", _du_dot_du = " <<_du_dot_du[_qp]<< std::endl;
  return _test[_i][_qp]*_phi[_j][_qp]*_du_dot_du[_qp];
}

_qp: 0, _du_dot_du = -0.154701
_qp: 1, _du_dot_du = 1
_qp: 2, _du_dot_du = 1
_qp: 3, _du_dot_du = 2.1547

[andrsd]

(1) This is my testing input file:

[Mesh]
  type = GeneratedMesh
  dim = 2
  xmin = 0
  xmax = 1
  ymin = 0
  ymax = 1
  nx = 10
  ny = 10
[]

[Functions]
  [./ffn]
    type = ParsedFunction
    value = 1
  [../]

  [./exact]
    type = ParsedFunction
    value = t
  [../]
[]

[Variables]
  [./u]
  [../]
[]

[Kernels]
  [./td]
    type = TimeDerivative
    variable = u
  [../]

  [./uff]
    type = UserForcingFunction
    variable = u
    function = ffn
  [../]
[]

[BCs]
  [./all]
    type = FunctionDirichletBC
    variable = u
    function = exact
    boundary = 'left right top bottom'
  [../]
[]

[Executioner]
  type = Transient
  dt = 0.1
  num_steps = 1

  solve_type = 'NEWTON'
[]

[Outputs]
  exodus = true
[]

So the only contribution to the Jacobian matrix is from the time derivative kernel. You can run it with -snes_type test is you want to, but I typically do that with -snes_check_jacobian -snes_check_jacobian_view.

(2) The residual printout by itself is not enough, I do agree, but since I am solving a linear problem, it has to solve like shown. Moreover, the important parts are the lines with 1.77973e-10 = ||J - Jfd||//J|| 1.12845e-09 = ||J - Jfd||. Those are the norms of the matrices, where J is what MOOSE computes and Jfd is the finite difference version. So you can see they are good.

(3) You have Diffusion and Reaction kernel on, so they contribute to the matrix as well, but I guess their Jacobians are exact, because they are simple kernels.

(4) The differences in your problem are:

• Constant monomials
• Adaptivity on

Those might be sources of possible problems.

I did not try to run your problem yet, but I will do it soon and will get back to you...

[lw4992]

I follow you test. It's OK both with and without "fix". The difference of our test is your test using Lagrange basis function while mine is monomials. If modifying your test by monomials, the issue would occur. Oops, I don't know why!

[lw4992]

I know why! Because for Lagrange basis, \sum phi[i] = 1 is always right, so _du_dot_du[_qp]=1/dt as MooseVariable.C, but for monomials, \sum phi[i] != 1, _du_dot_du[_qp] != 1/dt, it's wrong! Yes, it is the reason!

[YaqiWang]

These seems explaining things. We are lucky not having lots of transient tests using shape functions other than Lagrange. Question is that what are those 25 failing tests?

[andrsd]

@YaqiWang It does not matter what they were, the patch I used was wrong (as I said above).

Also, I know how to fix this problem. And it is not the way it was proposed (just for the record).

[YaqiWang]

@andrsd did you see any failing tests with the correct fix? Sorry for my blindness not noticing what you have said.

[andrsd]

@YaqiWang I am working on the correct fix ;-) The proposed one is not the right way to fix stuff and I did not try to run tests with it, so I do not know.