Mesh dependent convergence rate and Jacobian

[Nathan Miller]

I'm running into something strange for the convergence rate of a simple test problem. Using the newton solver I can get quadratic convergence for cubes of 2x2x2 and 4x4x4 elements generated by MOOSE. When I attempt to go to an 8x8x8 cube however the solver is converging, at best, sub-linearly.

When I ran the jacobian analyzer with randomICs for the 4x4x4 case I found that some of the values were significantly wrong. When running it for the 2x2x2 case I'm off by at most O 1e-2%.

Has anyone observed this sort of behavior or have any suggestions for how to debug?

[Alexander Lindsay]

Do you get good convergence if you use the FDP preconditioner with full = true? (Assuming that you're not using any constraints...)

--
You received this message because you are subscribed to the Google Groups "moose-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to moose-users+unsubscribe@googlegroups.com.
Visit this group at https://groups.google.com/group/moose-users.
To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/7f21c788-b0f4-44e9-b6c4-b9d71d2ca5da%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Nathan Miller]

I tried running with the finite difference preconditioner and I still don't see any convergence for an 8x8x8 problem.

On Monday, May 7, 2018 at 2:54:29 PM UTC-6, Alexander Lindsay wrote:

Do you get good convergence if you use the FDP preconditioner with full = true? (Assuming that you're not using any constraints...)

On Mon, May 7, 2018 at 2:21 PM, Nathan Miller <nathan...@gmail.com> wrote:

I'm running into something strange for the convergence rate of a simple test problem. Using the newton solver I can get quadratic convergence for cubes of 2x2x2 and 4x4x4 elements generated by MOOSE. When I attempt to go to an 8x8x8 cube however the solver is converging, at best, sub-linearly.

When I ran the jacobian analyzer with randomICs for the 4x4x4 case I found that some of the values were significantly wrong. When running it for the 2x2x2 case I'm off by at most O 1e-2%.

Has anyone observed this sort of behavior or have any suggestions for how to debug?

--
You received this message because you are subscribed to the Google Groups "moose-users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to moose-users...@googlegroups.com.

[Nathan Miller]

I'm not using any constraints either.

Have you ever seen something like this?

On Monday, May 7, 2018 at 2:54:29 PM UTC-6, Alexander Lindsay wrote:

Do you get good convergence if you use the FDP preconditioner with full = true? (Assuming that you're not using any constraints...)

On Mon, May 7, 2018 at 2:21 PM, Nathan Miller <nathan...@gmail.com> wrote:

I'm running into something strange for the convergence rate of a simple test problem. Using the newton solver I can get quadratic convergence for cubes of 2x2x2 and 4x4x4 elements generated by MOOSE. When I attempt to go to an 8x8x8 cube however the solver is converging, at best, sub-linearly.

When I ran the jacobian analyzer with randomICs for the 4x4x4 case I found that some of the values were significantly wrong. When running it for the 2x2x2 case I'm off by at most O 1e-2%.

Has anyone observed this sort of behavior or have any suggestions for how to debug?

--
You received this message because you are subscribed to the Google Groups "moose-users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to moose-users...@googlegroups.com.

[Derek Gaston]

Can you show us what your convergence history looks like?  I suspect that you're not getting convergence in your linear solver...

Derek

--

[Nathan Miller]

No problem. Here is the output from Moose. Everything seems to be converging in the linear solver but the non-linear rates aren't quadratic like expected. Note that this is for my hand-coded jacobian. I'll set the finite difference jacobian to run overnight just to verify.

Framework Information:
MOOSE Version:           git commit 1ad60ce on 2018-03-26
LibMesh Version:         727d0685dc96e841c5805608653410ea0cd926a6
PETSc Version:           3.7.6
Current Time:            Mon May  7 17:05:17 2018
Executable Timestamp:    Mon May  7 09:55:02 2018

Parallelism:
  Num Processors:          1
  Num Threads:             1

Mesh:
  Parallel Type:           replicated
  Mesh Dimension:          3
  Spatial Dimension:       3
  Nodes:                  
    Total:                 729
    Local:                 729
  Elems:                  
    Total:                 512
    Local:                 512
  Num Subdomains:          1
  Num Partitions:          1

Nonlinear System:
  Num DOFs:                8748
  Num Local DOFs:          8748
  Variables:               { "disp_x" "disp_y" "disp_z" "phi_xx" "phi_yy" "phi_zz" "phi_yz" "phi_xz"
                             "phi_xy" "phi_zy" "phi_zx" "phi_yx" }
  Finite Element Types:    "LAGRANGE"
  Approximation Orders:    "FIRST"

Execution Information:
  Executioner:             Steady
  Solver Mode:             NEWTON



 0 Nonlinear |R| = 2.019572e+05
      0 Linear |R| = 2.019572e+05
    0 KSP unpreconditioned resid norm 2.019572436327e+05 true resid norm 2.019572436327e+05 ||r(i)||/||b|| 1.000000000000e+00
      1 Linear |R| = 5.643791e-10
    1 KSP unpreconditioned resid norm 5.643790593312e-10 true resid norm 4.676469823410e-10 ||r(i)||/||b|| 2.315574197434e-15
Iteratively computed extreme singular values: max 1. min 1. max/min 1.
 1 Nonlinear |R| = 9.020063e+04
      0 Linear |R| = 9.020063e+04
    0 KSP unpreconditioned resid norm 9.020062976541e+04 true resid norm 9.020062976541e+04 ||r(i)||/||b|| 1.000000000000e+00
      1 Linear |R| = 2.193070e-10
    1 KSP unpreconditioned resid norm 2.193069613875e-10 true resid norm 2.582244187682e-10 ||r(i)||/||b|| 2.862778446667e-15
Iteratively computed extreme singular values: max 1. min 1. max/min 1.
 2 Nonlinear |R| = 3.721753e+04
      0 Linear |R| = 3.721753e+04
    0 KSP unpreconditioned resid norm 3.721753001906e+04 true resid norm 3.721753001906e+04 ||r(i)||/||b|| 1.000000000000e+00
      1 Linear |R| = 6.079513e-11
    1 KSP unpreconditioned resid norm 6.079512888811e-11 true resid norm 7.850599231781e-11 ||r(i)||/||b|| 2.109382118523e-15
Iteratively computed extreme singular values: max 1. min 1. max/min 1.
 3 Nonlinear |R| = 1.466025e+04
      0 Linear |R| = 1.466025e+04
    0 KSP unpreconditioned resid norm 1.466024951488e+04 true resid norm 1.466024951488e+04 ||r(i)||/||b|| 1.000000000000e+00
      1 Linear |R| = 1.291423e-10
    1 KSP unpreconditioned resid norm 1.291422657800e-10 true resid norm 1.538891688634e-10 ||r(i)||/||b|| 1.049703613211e-14
Iteratively computed extreme singular values: max 1. min 1. max/min 1.
 4 Nonlinear |R| = 1.334611e+04
      0 Linear |R| = 1.334611e+04
    0 KSP unpreconditioned resid norm 1.334611194484e+04 true resid norm 1.334611194484e+04 ||r(i)||/||b|| 1.000000000000e+00
      1 Linear |R| = 5.210927e-10
    1 KSP unpreconditioned resid norm 5.210927243260e-10 true resid norm 9.725450037745e-10 ||r(i)||/||b|| 7.287103598368e-14
Iteratively computed extreme singular values: max 1. min 1. max/min 1.
 5 Nonlinear |R| = 5.320058e+03
      0 Linear |R| = 5.320058e+03
    0 KSP unpreconditioned resid norm 5.320058468310e+03 true resid norm 5.320058468310e+03 ||r(i)||/||b|| 1.000000000000e+00
      1 Linear |R| = 2.562561e-10
    1 KSP unpreconditioned resid norm 2.562560752959e-10 true resid norm 2.962416132087e-10 ||r(i)||/||b|| 5.568390177914e-14
Iteratively computed extreme singular values: max 1. min 1. max/min 1.
 6 Nonlinear |R| = 4.955342e+03
      0 Linear |R| = 4.955342e+03
    0 KSP unpreconditioned resid norm 4.955342446258e+03 true resid norm 4.955342446258e+03 ||r(i)||/||b|| 1.000000000000e+00
      1 Linear |R| = 1.738264e-10
    1 KSP unpreconditioned resid norm 1.738264209378e-10 true resid norm 1.703872650622e-10 ||r(i)||/||b|| 3.438455907136e-14
Iteratively computed extreme singular values: max 1. min 1. max/min 1.
 7 Nonlinear |R| = 4.134356e+03
      0 Linear |R| = 4.134356e+03
    0 KSP unpreconditioned resid norm 4.134355794340e+03 true resid norm 4.134355794340e+03 ||r(i)||/||b|| 1.000000000000e+00
      1 Linear |R| = 1.541319e-10
    1 KSP unpreconditioned resid norm 1.541318668088e-10 true resid norm 2.064275850275e-10 ||r(i)||/||b|| 4.992980655175e-14
Iteratively computed extreme singular values: max 1. min 1. max/min 1.
 8 Nonlinear |R| = 3.805806e+03
      0 Linear |R| = 3.805806e+03
    0 KSP unpreconditioned resid norm 3.805806093956e+03 true resid norm 3.805806093956e+03 ||r(i)||/||b|| 1.000000000000e+00
      1 Linear |R| = 2.281632e-10
    1 KSP unpreconditioned resid norm 2.281632450583e-10 true resid norm 5.322363730831e-10 ||r(i)||/||b|| 1.398485261581e-13
Iteratively computed extreme singular values: max 1. min 1. max/min 1.
 9 Nonlinear |R| = 3.486940e+03

[Derek Gaston]

It looks like you're using a direct solver.  Unless this is extremely nonlinear... I would say the culprit is in your Jacobian.  If you switch to PJFNK does it fix it?

Derek

To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/ad991787-be29-4950-8889-fa3ebb921c07%40googlegroups.com.

[Nathan Miller]

I've tried both Newton and PJFNK and neither is able to converge in this case.

There does seem to be a problem with my Jacobian which only appears when I start using more elements. I'm using a Material to compute my stress measures. I am storing these material properties as std::vectors. Could the resizing of these vectors be a source of the problem somehow? I could replace these with fixed size arrays without too much difficulty.

[Nathan Miller]

I can see from the following thread that vector properties are okay: https://groups.google.com/forum/#!topic/moose-users/2Hz1YiEjnZ4

I communicate the gradients of my stress measures through vectors of vectors (std::vector<std::vector<double>>) which I use to construct the Jacobians of my kernels w.r.t. the degrees of freedom. I've declared my vector Material properties as double rather than Real. Is that an issue?

[Alexander Lindsay]

If FDP also fails then I think there's something more going on in this problem than potentially a poor hand-coded Jacobian.

Does it display poor convergence from the very beginning? 8000 DOFs is starting to get a little large for this...but you could try running with the petsc options (using NEWTON): `-pc_type svd -pc_svd_monitor`? An alternative would be to run with `-pc_type none -ksp_monitor_singular_value` (using either NEWTON or PJFNK). These options would tell us what the condition number of your system is. Especially if you're using a direct solver like Derek mentioned, you may be hiding issues in your linear solve.

How many linear iterations does it take to converge when using PJFNK? I see that you're using `-ksp_monitor_true_residual`. I would be curious to see the output of a linear solve with PJFNK.

Does this work involve custom kernels, bcs, etc? We might be able to help if you linked us to a branch (preferred) or attached an input file if we don't need any custom code.

Alex

To unsubscribe from this group and stop receiving emails from it, send an email to moose-users+unsubscribe@googlegroups.com.

Visit this group at https://groups.google.com/group/moose-users.

To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/30022614-f639-4aa9-bdc8-09d874d45fc3%40googlegroups.com.

[Nathan Miller]

FDP also failed using the Newton solver. Interestingly the first few nonlinear iterations have the same residual as with my hand-coded jacobian.

The convergence for 8x8x8 is bad from the very beginning while for smaller problems it is quadratic.

When I ran with -pc_type svd and -pc_svd_monitor using my hand-coded jacobian it stalled (as you indicated it might). My workstation is a bit old (and my jacobian expensive) so that is the likely issue. I'll keep it running in the background and update later.

With PJFNK I still find that the linear system doesn't converge on these larger problems. On a 2x2x2 grid using the following PETSc options:
[Preconditioning]
  [./SMP]
    type = SMP
    #type = FDP
    full = true
  [../]
[]

[Executioner]
  type = Steady
  solve_type = 'PJFNK'
  petsc_options = '-snes_ksp_ew -ksp_monitor_true_residual -ksp_compute_singularvalues'
  petsc_options_iname = '-pc_type -sub_pc_type -pc_asm_overlap -ksp_gmres_restart'
  petsc_options_value = 'asm      lu           1               101               '
[]

I get the following (2x2x2 grid):

 0 Nonlinear |R| = 3.657852e+03
      0 Linear |R| = 3.657852e+03
    0 KSP unpreconditioned resid norm 3.657851852883e+03 true resid norm 3.657851852883e+03 ||r(i)||/||b|| 1.000000000000e+00
      1 Linear |R| = 1.066267e-04
    1 KSP unpreconditioned resid norm 1.066266618933e-04 true resid norm 1.066266618196e-04 ||r(i)||/||b|| 2.915007663187e-08

Iteratively computed extreme singular values: max 1. min 1. max/min 1.

 1 Nonlinear |R| = 2.198872e+02
      0 Linear |R| = 2.198872e+02
    0 KSP unpreconditioned resid norm 2.198871848142e+02 true resid norm 2.198871848142e+02 ||r(i)||/||b|| 1.000000000000e+00
      1 Linear |R| = 2.207833e-05
    1 KSP unpreconditioned resid norm 2.207833388103e-05 true resid norm 2.207833388140e-05 ||r(i)||/||b|| 1.004075517182e-07

Iteratively computed extreme singular values: max 1. min 1. max/min 1.

 2 Nonlinear |R| = 2.842126e+00
      0 Linear |R| = 2.842126e+00
    0 KSP unpreconditioned resid norm 2.842126011735e+00 true resid norm 2.842126011735e+00 ||r(i)||/||b|| 1.000000000000e+00
      1 Linear |R| = 7.695418e-07
    1 KSP unpreconditioned resid norm 7.695417882062e-07 true resid norm 7.695417881913e-07 ||r(i)||/||b|| 2.707627265694e-07

Iteratively computed extreme singular values: max 1. min 1. max/min 1.

 3 Nonlinear |R| = 1.843390e-03
      0 Linear |R| = 1.843390e-03
    0 KSP unpreconditioned resid norm 1.843390298813e-03 true resid norm 1.843390298813e-03 ||r(i)||/||b|| 1.000000000000e+00
      1 Linear |R| = 7.104829e-10
    1 KSP unpreconditioned resid norm 7.104828873693e-10 true resid norm 7.104828873698e-10 ||r(i)||/||b|| 3.854218435604e-07

Iteratively computed extreme singular values: max 1. min 1. max/min 1.

 4 Nonlinear |R| = 1.947349e-09

But on a 8x8x8 grid I get:

 0 Nonlinear |R| = 2.019572e+05
      0 Linear |R| = 2.019572e+05
    0 KSP unpreconditioned resid norm 2.019572436327e+05 true resid norm 2.019572436327e+05 ||r(i)||/||b|| 1.000000000000e+00

      1 Linear |R| = 5.441531e-02
    1 KSP unpreconditioned resid norm 5.441531418677e-02 true resid norm 5.441531418322e-02 ||r(i)||/||b|| 2.694397745008e-07

Iteratively computed extreme singular values: max 1. min 1. max/min 1.

 1 Nonlinear |R| = 9.020061e+04
      0 Linear |R| = 9.020061e+04
    0 KSP unpreconditioned resid norm 9.020060927729e+04 true resid norm 9.020060927729e+04 ||r(i)||/||b|| 1.000000000000e+00
      1 Linear |R| = 4.024374e-03
    1 KSP unpreconditioned resid norm 4.024373994073e-03 true resid norm 4.024373994732e-03 ||r(i)||/||b|| 4.461581830740e-08

Iteratively computed extreme singular values: max 1. min 1. max/min 1.

 2 Nonlinear |R| = 3.721752e+04
      0 Linear |R| = 3.721752e+04
    0 KSP unpreconditioned resid norm 3.721752052275e+04 true resid norm 3.721752052275e+04 ||r(i)||/||b|| 1.000000000000e+00
      1 Linear |R| = 8.978024e-04
    1 KSP unpreconditioned resid norm 8.978023858695e-04 true resid norm 8.978023859251e-04 ||r(i)||/||b|| 2.412311119372e-08

Iteratively computed extreme singular values: max 1. min 1. max/min 1.
 3 Nonlinear |R| = 1.466025e+04

So as you can see it is following the Newton nonlinear residual exactly though it still doesn't converge. When I use boomerAMG it isn't able to solve the linear problem at all.

This branch does involve custom kernels which call a shared object library compiled external to MOOSE. Let me do some checking about sharing my repositories. I may have some restrictions there.

Nathan

[Daniel Schwen]

> I communicate the gradients of my stress measures through vectors of vectors (std::vector<std::vector<double>>) which I use to construct the Jacobians of my kernels

Can we actually see the code? Otherwise this is just poking in the dark.

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/6a1aef62-0124-4fc3-84ef-45a44ce05d18%40googlegroups.com.

[Peterson, JW]

On Tue, May 8, 2018 at 11:13 AM, Nathan Miller <nathan...@gmail.com> wrote:

FDP also failed using the Newton solver. Interestingly the first few nonlinear iterations have the same residual as with my hand-coded jacobian.

Nathan, is it possible for you to push your code somewhere that we can take a look at it? Particularly the Material properties you've written, but also the Kernels.

--

John

[Nathan Miller]

Sure, I just needed to check that I was cleared to do that.

You will need two repositories. One for my MOOSE app and the other for the shared libraries my App calls. The ReadMe in both repositories should show how to set it up.

The repository for my MOOSE app (tardigrade) is located here: https://bitbucket.org/NateAM2/tardigrade/src/master/
The repository for my shared libaries is located here:                  https://bitbucket.org/NateAM2/micromorphic_element/src/master/

You will also need the Eigen library if you want to try and compile it (http://eigen.tuxfamily.org/).

All of the runs seen above were generated with: tardigrade/test/tests/kernels/stretch_y
Note there are no BCs on my dof tensor phi. You shouldn't need that because of the coupling with disp.

Thanks for your help,

Nathan

[Peterson, JW]

On Tue, May 8, 2018 at 11:46 AM, Nathan Miller <nathan...@gmail.com> wrote:

Sure, I just needed to check that I was cleared to do that.

You will need two repositories. One for my MOOSE app and the other for the shared libraries my App calls. The ReadMe in both repositories should show how to set it up.

The repository for my MOOSE app (tardigrade) is located here: https://bitbucket.org/NateAM2/tardigrade/src/master/
The repository for my shared libaries is located here:                  https://bitbucket.org/NateAM2/micromorphic_element/src/master/

You will also need the Eigen library if you want to try and compile it (http://eigen.tuxfamily.org/).

All of the runs seen above were generated with: tardigrade/test/tests/kernels/stretch_y
Note there are no BCs on my dof tensor phi. You shouldn't need that because of the coupling with disp.

I'm taking a quick look now. One thing that would make your input file a lot shorter is if you add a [GlobalParams] section with all the coupling, i.e.

[GlobalParams]

    u1     = disp_x

    u2     = disp_y

    u3     = disp_z

    phi_11 = phi_xx

    phi_22 = phi_yy

    phi_33 = phi_zz

    phi_23 = phi_yz

    phi_13 = phi_xz

    phi_12 = phi_xy

    phi_32 = phi_zy

    phi_31 = phi_zx

    phi_21 = phi_yx

[]

then you don't have to repeat it in every Kernel.

--

John

[Nathan Miller]

Great! Thanks!

I will add the [GlobalParams] block. That definitely helps a lot.

Alexander requested earlier that I do a run with -pc_type svd -pc_svd_monitor. The results of the first iteration are:

 0 Nonlinear |R| = 2.019572e+05

      SVD: condition number 4.069453376632e+09, 0 of 8748 singular values are (nearly) zero
      SVD: smallest singular values: 1.118221874211e-03 1.235165226648e-03 1.267933471671e-03 1.302050300766e-03 1.615674918363e-03
      SVD: largest singular values : 3.823772593004e+06 3.994559437820e+06 4.263501040694e+06 4.263501040695e+06 4.550551781832e+06

      0 Linear |R| = 2.019572e+05
    0 KSP unpreconditioned resid norm 2.019572436327e+05 true resid norm 2.019572436327e+05 ||r(i)||/||b|| 1.000000000000e+00

      1 Linear |R| = 1.373937e-08
    1 KSP unpreconditioned resid norm 1.373936624825e-08 true resid norm 1.406425777771e-08 ||r(i)||/||b|| 6.963977882014e-14

Iteratively computed extreme singular values: max 1. min 1. max/min 1.
 1 Nonlinear |R| = 9.020063e+04

[Daniel Schwen]

> SVD: condition number 4.069453376632e+09

Uhm, I'm not a mathematician, but that condition number strikes my as a little high...

--

You received this message because you are subscribed to the Google Groups "moose-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to moose-users...@googlegroups.com.
Visit this group at https://groups.google.com/group/moose-users.

To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/21c8cd08-0642-4fd6-9f06-6b0104ad5f32%40googlegroups.com.

[Peterson, JW]

On Tue, May 8, 2018 at 2:03 PM, Daniel Schwen <dan...@schwen.de> wrote:

> SVD: condition number 4.069453376632e+09

Uhm, I'm not a mathematician, but that condition number strikes my as a little high...

This condition number does seem high for a small test problem, but it doesn't guarantee that anything in particular is wrong, like a condition number of ~1e16 would.

That said, I'm also not sure that I buy the idea that you don't need *any* boundary conditions on the phi_xx, phi_yy, etc. equations. I don't really understand the mathematical character of the tensor equations, but they are PDEs, and those usually require BCs.

After looking through the code for a bit, I don't see anything terribly wrong as far as MOOSE coding practices. (I am a bit surprised that the getFunction() calls in MicromorphicMaterial don't complain when you don't provide the Functions in question, but that might be OK as long as you don't actually try to call them...)

Most of the complexity of the app is obviously buried in the balance_equations:: namespace, which I did not look through in detail. But it seems to me that when calling interfaces like this:

balance_equations::compute_internal_couple_jacobian(_component_i,         _component_j,      _dof_num, 

                                                        _test[_i][_qp],        dNdx,             _phi[_j][_qp],          detadx,

                                                        _DcauchyDgrad_u[_qp], _DcauchyDphi[_qp], _DcauchyDgrad_phi[_qp],

                                                        _DsDgrad_u[_qp],      _DsDphi[_qp],      _DsDgrad_phi[_qp],

                                                        _DmDgrad_u[_qp],      _DmDphi[_qp],      _DmDgrad_phi[_qp],

                                                        dcdUint);

it would be easy to get one or two things in the wrong order... if you are looking for a place to start debugging, that might be a good one.

I also have not seen what you are using for the weak form of your equations (sorry if already posted) but it would be good if you could point us to that as well.

--

John

[Nathan Miller]

The condition number does seem to be a problem.

I met with Jed Brown today (CU professor and PETSc developer) and he helped me break the problem down a bit further. It looks like some of my diagonal terms in my Jacobian are significantly different than what I would have expected compared to the other terms (1e4 vs 1e0). So I think there may be a problem there. I'll dig into that this evening.

In terms of the getFunction() calls, is there a better way to pull in those functions? I'm using that for my method of manufactured solutions verification.

In terms of the balance equations, the weak form is given in the figure below. I like to use index notation so that is what I generally use. Note that in this context A is the cauchy stress (sigma), B is the symmetric stress (s) and C is the higher order stress (m) a_i are my displacements and b_ij are my tensoral phi values. If you are looking for a write-up which gives a far more detailed derivation that can be found here: http://ceae.colorado.edu/~regueiro/publications/reports/Regueiro_ARO_report_final_08080.pdf

[Nathan Miller]

Correction on my last line. I edited my graphics to reflect the notation used in the paper but didn't update my text before posting.

To recap: "sigma" is the cauchy stress, "s" is the symmetric stress, and "m" is the higher order stress. The other terms are not implemented into kernels yet though they exist in balance_equations.cpp.

I'm still working through my Jacobian to try and understand where things are going wrong. I have some negative values which don't really make sense given that this problem shouldn't have instabilities. I'll check this against the finite difference version tomorrow. I might be mistaken but I was sure that I should be able to let my phi DOF float everywhere and the solution will resolve. Maybe I was incorrect though.

On a related note, you can see from my PDEs that one of them has units of force and the other has units of stress. In looking at my Jacobian, I found some of my diagonal terms would vary by as much as 3 orders of magnitude. I scaled each of my variables by the inverse average order of magnitude of each diagonal term of my Jacobian which made my diagonal terms order 1 but didn't really seem to help my ill-conditioned matrix. It could be that my off-diagonal terms were still poor. I'm still not sure that I don't have a bug somewhere but given that the finite difference and my hand-coded jacobian result in the same convergence profile I'm wondering if it has more to do with poor conditioning than anything else.

Are there any good resources for how to deal with ill-conditioned matrices?

[Alexander Lindsay]

I'm hoping I can try running your problem later today, but can't promise. In the meantime...how many processes are you running with? The fact that you're converging in one linear iteration makes me think that you're only running with one process so only `-sub_pc_type` is being used... You've said that boomerAMG doesn't work. Can you try with `-pc_type ilu`? Also if you're running with `-ksp_monitor_true_residual` you might as well specify `print_linear_residuals = false` in your Outputs block to remove some of that output clutter.

On your 2x2x2 and 4x4x4 cases that have solved, have you looked at the outputs? Do they look like they should? It sounds like you have some MMS in your code...have you verified against MMS? This knowledge would help us address concern about missing boundary conditions.

--
You received this message because you are subscribed to the Google Groups "moose-users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to moose-users+unsubscribe@googlegroups.com.

Visit this group at https://groups.google.com/group/moose-users.

To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/a1c43ade-0f70-41b7-b27b-665b2db3e7b4%40googlegroups.com.

[Nathan Miller]

That would be great! I think you will be the first to try and do that so I'll look forward to hearing your feedback. I've pushed some updates for the stretch_y.i file so make sure you pull a fresh copy. In particular, some of the material parameters I was using were sub-optimal. In particular, they were contributing to the ill-conditioned nature of my jacobian. I've defined the internal force from the balance of linear momentum as - _grad_test[_i][_qp] * cauchy[_qp] which has caused some of my diagonal entries to my jacobian to be negative.

I am only running with one process currently. I wanted to try and limit potential sources of error. I tried using -pc_type ilu and -print_linear_residuals false using the finite difference preconditioner to eliminate that source of error. I've attached the results of that if you want to see them. The nonlinear convergence was linear in this case for the 4x4x4 problem and sublinear (though it converged!) for the 8x8x8 problem. The ratio of my max and min singular values got to 3184.17 by the end so it was starting to get a bit ill-conditioned. I would be happy for any suggestions for how to apply scaling to try and fix the ill-conditioning and get this to converge in fewer iterations. Using asm instead of ilu didn't get convergence with the 8x8x8 case.

My MMS solutions do give me what I expect though I need to add more to exercise all of the parts of the model that could be giving problems. Right now they are just constant shifts in u_i and phi_ij so everything should be real boring. I'll add linear shifts next to see what happens.

In speaking to my advisor, we should be able to allow phi to float without there being any issues for the linear-elastic case. Also, there is a sign error in the strong form I posted above. I'm so used to dealing with the weak form.

[Peterson, JW]

On Wed, May 9, 2018 at 10:58 AM, Nathan Miller <nathan...@gmail.com> wrote:

That would be great! I think you will be the first to try and do that so I'll look forward to hearing your feedback. I've pushed some updates for the stretch_y.i file so make sure you pull a fresh copy. In particular, some of the material parameters I was using were sub-optimal. In particular, they were contributing to the ill-conditioned nature of my jacobian. I've defined the internal force from the balance of linear momentum as - _grad_test[_i][_qp] * cauchy[_qp] which has caused some of my diagonal entries to my jacobian to be negative.

I am only running with one process currently. I wanted to try and limit potential sources of error. I tried using -pc_type ilu and -print_linear_residuals false using the finite difference preconditioner to eliminate that source of error. I've attached the results of that if you want to see them. The nonlinear convergence was linear in this case for the 4x4x4 problem and sublinear (though it converged!) for the 8x8x8 problem. The ratio of my max and min singular values got to 3184.17 by the end so it was starting to get a bit ill-conditioned. I would be happy for any suggestions for how to apply scaling to try and fix the ill-conditioning and get this to converge in fewer iterations. Using asm instead of ilu didn't get convergence with the 8x8x8 case.

I wouldn't say O(10^3) is particularly ill-conditioned. Note that the condition number grows like h^{-2} for second-order problems, so the Jacobian will always become more ill-conditioned as you refine the mesh.

The usual approaches for fixing ill-conditioning are equation scaling (easy, but may not completely fix the problem) and developing a consistent non-dimensionalized formulation of your problem (typically involves more work). 

Since the unknowns in your problem have different units/scales (forces vs. stresses) it might be worthwhile thinking about non-dimensionalization. This is more common in fluid mechanics applications, i.e. Reynolds number and other non-dimensional parameters will appear in the governing equations, I'm not sure how much it is used in tensor mechanics applications.

--

John

[Nathan Miller]

John,

I had been considering developing the non-dimensionalized form though I was hoping to be able to avoid that. I haven't heard of people doing it for solid mechanics either though I would expect someone has tried it.

Nathan

[Alexander Lindsay]

You said the 8x8x8 was sublinear, but the solve information you attached makes it look more like super-linear to me. Here's a plot of mu_k:

Maybe I'm misunderstanding

--

You received this message because you are subscribed to the Google Groups "moose-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to moose-users+unsubscribe@googlegroups.com.
Visit this group at https://groups.google.com/group/moose-users.

To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/f1f49c8d-63cb-44a9-87b8-fb4863cbbc1c%40googlegroups.com.

[Alexander Lindsay]

Moreover, using `-snes_ksp_ew` makes such convergence plots a little dubious to me since you're solving the linearized system by a different amount each iteration...

[Nathan Miller]

The attached info was for the 4x4x4 case not for the 8x8x8. Sorry for the confusion. I've attached both the 4x4x4 and the 8x8x8 plots.

When I don't use -snes_ksp_ew for the 4x4x4 case I get quadratic convergence using the article you linked for the second to last iteration. For the 8x8x8 case using ilu and without -snes_ksp_ew the solution starts out linear/sublinear but then starts doing much better towards the end. So maybe my question really is how do I get quadratic convergence the whole way as opposed to just during the final few iterations.

To unsubscribe from this group and stop receiving emails from it, send an email to moose-users...@googlegroups.com.

Visit this group at https://groups.google.com/group/moose-users.

[Peterson, JW]

On Wed, May 9, 2018 at 12:51 PM, Nathan Miller <nathan...@gmail.com> wrote:

The attached info was for the 4x4x4 case not for the 8x8x8. Sorry for the confusion. I've attached both the 4x4x4 and the 8x8x8 plots.

When I don't use -snes_ksp_ew for the 4x4x4 case I get quadratic convergence using the article you linked for the second to last iteration. For the 8x8x8 case using ilu and without -snes_ksp_ew the solution starts out linear/sublinear but then starts doing much better towards the end. So maybe my question really is how do I get quadratic convergence the whole way as opposed to just during the final few iterations.

Quadratic convergence in Newton is really only guaranteed when you are near the root...

--

John

[Nathan Miller]

Of course. I was just surprised to see such a difference between the two mesh densities. Maybe that is more typical than I thought.

[Derek Gaston]

On Wed, May 9, 2018 at 12:51 PM Nathan Miller <nathan...@gmail.com> wrote:

So maybe my question really is how do I get quadratic convergence the whole way as opposed to just during the final few iterations.

 Several options:

1.  Take a better initial guess (especially try to get your initial value for your BCs to match - if you are using DirichletBCs switch to PresetBC instead).

2.  "ramp up" a parameter over "time".  Switch to Transient and increase your load every timestep (use a time dependent BC or forcing function).

3.  Pseudo-transient simulation to steady state.  Add a time derivative term and turn on the steady state checking for Transient... and iterate to steady state.  As a bonus use "predictor_scale" to take even better guesses each timestep.

4.  Do #4 but use Explicit to start the simulation (using lump-preconditioned ActuallyExplicitEuler)... then switch to implicit with an adaptive timestepper to keep taking bigger steps as you get close to the steady state solution.  Could do this with restart (although we are planning on making this possible in general).

With Newton... it's all about solving from where you are to the root.  If the root you are looking for is closer to where you're starting from - then you will get better solver performance.

By using the above techniques you're making it so that the Newton solver doesn't have very "far" to go... which should make it easier to get quadratic convergence.

Derek

[Nathan Miller]

Ah, okay. I didn't consider using a pseudo time. I will give that a shot.

[Derek Gaston]

I would recommend using CoefTimeDerivative so that you can balance the strength of the time-derivative vs everything else.  A "stronger" time-derivative will mean that things will change more slowly... which will make the solve easy... but if you go too far then you are wasting cpu cycles doing a bunch of unnecessary solves.

The equivalent thing is to change your timestep size.  Smaller timesteps will be easier to solve... with the same tradeoff.

I find that when I'm doing pseudo transients I typically like to take timesteps around ~1... so I adjust CoefTimeDerivative to give me a good amount of movement with a dt of ~1... then I can let an adaptive timestepper (like SolutionTimeAdaptive or IterationAdaptiveDT) take over from there to get to the steady state as fast as possible.

Derek

--

You received this message because you are subscribed to the Google Groups "moose-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to moose-users...@googlegroups.com.
Visit this group at https://groups.google.com/group/moose-users.

To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/cff95404-e323-4ddb-b29b-4c9ca40d995a%40googlegroups.com.

[Derek Gaston]

Oooh - also... of course... you can adjust (or leave out!) the CoefTimeDerivative _per equation_.  You might have some "easy" equations that can take a "larger" change... then let the harder physics evolve more slowly.

All up to you.

Derek

[Nathan Miller]

I should have thought of using the transient solver to make the problem easier. That is exactly what I was doing in Abaqus without even thinking about it. So much is handled with that code behind the scenes that it is easy (for me) to forget exactly what is going on.

The addition of the transient solve fixed all of my problems with convergence.

Thanks to everyone for their help and patience. It really is appreciated.

Nathan

[Derek Gaston]

Great!  Thanks for letting us know.

Derek

To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/b0decea4-93e9-4600-939b-613567e52ed9%40googlegroups.com.