Solver setting for Poisson-Nernst-Planck (PNP) equations for electrochemical systems

714 views
Skip to first unread message

Shohei Ogawa

unread,
May 18, 2018, 4:25:22 PM5/18/18
to moose-users

Dear Moose developers,


I am working on the modeling of a fuel cell electrode in 3D by using Poisson-Nernst-Planck equations. So far I am studying for steady-state.

I have the Nernst-Planck equation is for ion concentrations including the diffusion and electromigration effect:

where c_i and Psi are the ion concentration and the electric potential in an electrolyte, respectively. I have a boundary condition for the ion species for chemical reaction, which exponentially depends on the electric potential and has some linear dependence on the ion concentration.


The electric potential appears in the Poisson's equation which is: 


The problem may be challenging because

1) The ion movement depends on the gradient of the electric potential

2) The two equations are coupled to each other.

3) The boundary condition is highly non-linear.


Has anyone worked on a similar model, and 


I have made sure my Jacobian expressions are correct for the kernel for the weak form of nabla * (F*z_i*mu_i*c_i nabla *Psi) and the boundary condition expressions. 

BoomerAMG in Hypre didn't decrease the linear residual at all. ASM worked better but takes thousands of iterations and the non-linear residual decrease is very small.


I am looking through the PETSc document (http://www.mcs.anl.gov/petsc/petsc-current/docs/), but also I would appreciate it if you could give me any suggestion!


Thank you,

Shohei


Derek Gaston

unread,
May 18, 2018, 4:55:30 PM5/18/18
to moose...@googlegroups.com
I would expect AMG to work pretty well here - depending  on the strength of that convection term.

1.  Try to solve a very small problem first
2.  Make sure your Jacobian is correct.
3.  Use LU to be sure that if you perfectly invert the Jacobian that everything works correctly
4.  Use FDP with LU to be doubly sure.

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/aac527c0-9399-43a0-a637-149494dfd293%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Derek Gaston

unread,
May 18, 2018, 4:57:17 PM5/18/18
to moose...@googlegroups.com

Shohei Ogawa

unread,
May 18, 2018, 5:10:29 PM5/18/18
to moose...@googlegroups.com
Hi Derek,
 
 Thank you very much for a quick reply and instructions.
 I will work on a simple problem first. If I can solve the problem with AMG, that would be great.

Regarding with the instructions,
3.  Use LU to be sure that if you perfectly invert the Jacobian that everything works correctly
4.  Use FDP with LU to be doubly sure.

Could you please give me more instructions? For example, what PETSc options should I have in an executioner block and what should I look at to make sure everything is correct?

Thank you,
Shohei


Anil Kunwar

unread,
May 18, 2018, 9:48:52 PM5/18/18
to moose-users
Dear Shohei,
I have developed the Nernst-Planck Electromigration term for the first PDE.
You can  have an overview of it at:

One of the things you may need to do in multiphysics coupling is the scaling of variables c_{i} and phi_{i}.
[Variables]
[./variable-1]
order = FIRST
family = LAGRANGE
initial_condition = 1.0
scaling = 1.0e+02
[../]
[./variable-2]
order = FIRST
family = LAGRANGE
initial_condition = 2.0
scaling = 1.0e+05
[../]
[]
The scaling of these variables must be done in such a way that the maximum residual should be in the same or 1 order different.

Yours Sincerely,
Anil Kunwar

Derek Gaston

unread,
May 18, 2018, 10:27:35 PM5/18/18
to moose...@googlegroups.com
Thanks Anil!  Always helpful to get physics experts chiming in!

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/be730973-b567-4a72-8e16-06fec56ebb03%40googlegroups.com.

Shohei Ogawa

unread,
May 19, 2018, 4:46:27 PM5/19/18
to moose-users
Anil,

 Thank you for introducing me your work! What preconditioner did you use to solve a problem including that kernel? I didn't find any input file in your repository which includes the kernel.

 Regarding with the scaling, I agree. The order of the residuals can be very different when we have multiple ions in particular.

Thank you,
Shohei 

Anil Kunwar

unread,
May 19, 2018, 8:27:12 PM5/19/18
to moose-users
Hi Shohei,
1. I had developed the NernstPlanckConvection kernel for analysing electrochemical migration examples, and it is still incomplete, as there is no any input file electrochemicalmigration directory.

2. It is possible to write an input file by yourself from a close matching example of electromigration (https://github.com/anilkunwar/danphe/blob/working_codes/examples/Electromigration/cu_transport_j_l1243_da.i), where the ConstantTensorElectricConvection utilizes constant current density vector. When you use NernstPlanckConvection kernel in place of ConstantTensorElectricConvection, then the voltage variable (secondary variable or coupled variable  as ion concentration is the primary variable 'u' in MOOSE terminology ) , must be introduced with additional kernels (kernels reresenting the second equation of your first post).  To see a similar example , with c variable (c) and temperature variable (T), please find https://github.com/anilkunwar/danphe/blob/working_codes/examples/Thermomigration/cu_transport_T250_l200.i . In this equation, T variable is coupled  to the ThermalConvection kernel by first introducing in the MatDiffusion kernel (with T variable). There is another MatDIffusion kernel (with c variable ) in the same input file.

3. And, notably, i guess, you first can try by developing simpler models ( with \nabla (\epsilon_r \epsilon_0) . \nabla (\phi) = 0). After it is solved, you can add the remaining two kernels of the Eq. 2 of your first post.

4. If you use the equation, with \nabla (\epsilon_r \epsilon_0) . \nabla (\phi) = 0, to represent the voltage distribution , in context of simpler model, the preconditioners that have been used in the above examples may work for you too. Now, if you add the time-derivate term (\partial \phi/ \partial t) and source term (term containing ion concentrations), the numerical implementation becomes more complicated ; and in that context , depending upon the difficulty faced, you may have to make extra workout on the preconditioner and executioner blocks.

Yours Sincerely,
Anil Kunwar

Anil Kunwar

unread,
May 19, 2018, 9:21:37 PM5/19/18
to moose-users
Hi Shohei,
Today, I have developed an example input file using NernstPlanckConvection kernel , and have just now uploaded in github repository, and it can be found in https://github.com/anilkunwar/danphe/blob/working_codes/examples/Electrochemicalmigration/cu_transport_phi_l450_small-length.i .
It models Cu species transport across Sn medium under the presence of steady state voltage difference.
This is a simple model with voltage distribution defined by \nabla. \sigma \nabla \phi = 0.
You can further extend or develop it by adding the time derivative and source term to the kernels of voltage variables. In source term of voltage variable, you have to couple concentration (as coupled variable).

Yours Sincerely,
Anil Kunwar

Shohei Ogawa

unread,
May 20, 2018, 3:38:42 PM5/20/18
to moose...@googlegroups.com
Hi Anil,

 Thank you very much for the detailed examples. I looked at your input files. Unfortunately, I wasn't able to `make` your example because of some errors. I believe other files in your local have to be pushed to the Github repository.

 I solved coupled Poisson-Nernst-Planck equation for a simple and already-solved problem and the result of my app seems to be correct. For this problem, Boomer AMG as a preconditioner was successful. I am wondering when I add my reaction term as a boundary condition, the convergence gets much more challenging. I will go back to simple model this week and try to see what is causing a problem.

Thank you,
Shohei

Shohei Ogawa

To unsubscribe from this group and stop receiving emails from it, send an email to moose-users+unsubscribe@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/b8839399-0c51-42ed-98d3-69e63263137f%40googlegroups.com.

anil kunwar

unread,
May 20, 2018, 9:50:51 PM5/20/18
to moose...@googlegroups.com
Dear Shohei,

-If you have compiled danphe app (danphe-working_codes branch), it should compile fine , and you must have to update some codes settings according to MOOSE version changes.
- I am happy to know that you are able to run a model in your own app.

Yours Sincerely,
Anil Kunwar

On Monday, May 21, 2018, 3:38:43 AM GMT+8, Shohei Ogawa <ogawa...@gmail.com> wrote:


Hi Anil,

 Thank you very much for the detailed examples. I looked at your input files. Unfortunately, I wasn't able to `make` your example because of some errors. I believe other files in your local have to be pushed to the Github repository.

 I solved coupled Poisson-Nernst-Planck equation for a simple and already-solved problem and the result of my app seems to be correct. For this problem, Boomer AMG as a preconditioner was successful. I am wondering when I add my reaction term as a boundary condition, the convergence gets much more challenging. I will go back to simple model this week and try to see what is causing a problem.

Thank you,
Shohei

Shohei Ogawa

On Sat, May 19, 2018 at 9:21 PM, Anil Kunwar <romagu...@gmail.com> wrote:
Hi Shohei,
Today, I have developed an example input file using NernstPlanckConvection kernel , and have just now uploaded in github repository, and it can be found in https://github.com/anilkunwar/ danphe/blob/working_codes/ examples/ Electrochemicalmigration/cu_ transport_phi_l450_small- length.i .

I am looking through the PETSc document (http://www.mcs.anl.gov/petsc/ petsc-current/docs/), but also I would appreciate it if you could give me any suggestion!


Thank you,

Shohei


--
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.co m.
Visit this group at https://groups.google.com/grou p/moose-users.
To view this discussion on the web visit https://groups.google.com/d/ms gid/moose-users/be730973-b567- 4a72-8e16-06fec56ebb03%40googl egroups.com.
For more options, visit https://groups.google.com/d/op tout.

--
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.

For more options, visit https://groups.google.com/d/ optout.

--
You received this message because you are subscribed to a topic in the Google Groups "moose-users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/moose-users/O5dLCn-oowE/unsubscribe.
To unsubscribe from this group and all its topics, send an email to moose-users...@googlegroups.com.

walkand...@gmail.com

unread,
May 21, 2018, 11:53:50 AM5/21/18
to moose-users
Thank you sincerely Anil,
this scaling is great helpful!

Thank you

在 2018年5月19日星期六 UTC+2上午3:48:52，Anil Kunwar写道：

Shohei Ogawa

unread,
May 21, 2018, 4:54:03 PM5/21/18
to moose-users
Hi Anil,

 Thank you. I just want to let you know that I was able to modify the programs a little and compile the app.

Shohei

Shohei Ogawa

Shohei Ogawa

unread,
May 21, 2018, 11:24:54 PM5/21/18
to moose-users
Hi Derek,

Following what you mentioned in the previous post below, I did some studies.

> 1.  Try to solve a very small problem first
I created a simple problem in 2D and solved with the ASM preconditioner. Basically, the result matched with the model implemented in COMSOL. In the future, with Moose, I want to solve the problem for a much larger domain and dense mesh in 3D for complex geometry.

> 2.  Make sure your Jacobian is correct.
I used the Jacobian debugging script and found no error so far.

> 3.  Use LU to be sure that if you perfectly invert the Jacobian that everything works correctly

petsc_options_value = 'lu'

petsc_options_iname = '-pc_type'

solve_type = PJFNK


Those setting in the executioner block worked. The simulation converged within a few non-linear iterations and a few linear iterations for each non-linear iteration.


> 4.  Use FDP with LU to be doubly sure.

petsc_options_value = 'lu'

petsc_options_iname = '-pc_type'

solve_type = FD


Those setting in the executioner block worked as well.


However, with

petsc_options_value = 'hypre boomeramg' petsc_options_iname = '-pc_type -pc_hypre_type' solve_type = PJFNK

The simulation stopped due to Linear solve did not converge due to DIVERGED_DTOL iterations 150

Much larger -ksp_gmres_restart value didn't work because the residual didn't decrease.


Shohei

Alexander Lindsay

unread,
May 22, 2018, 10:01:02 AM5/22/18
to moose...@googlegroups.com
If the convection term is strong, then multi-grid likely won't be effective

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/285923a4-ed05-4a94-bd6d-b56a36dc28d8%40googlegroups.com.

Kong, Fande

unread,
May 22, 2018, 10:36:08 AM5/22/18
to moose...@googlegroups.com
Just for clarification, AMG does not like non-elliptic equations much.  However, GMG may work for a wide range of applications including non-elliptic equations. We may be going to support this algorithm in the future.


Fande,

To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/CANFcJrFsPNgt_wdDTtD%2BQJ6SuAP5WOcoQQVKv3aMoCVDedHWgg%40mail.gmail.com.

Shohei Ogawa

unread,
May 22, 2018, 10:46:15 AM5/22/18
to moose...@googlegroups.com
Thank you for the comments! 
So is ASM the first choice now? Is there any good preconditioner which can work for a convection problem in parallel?

Shohei

--
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/CANFcJrFsPNgt_wdDTtD%2BQJ6SuAP5WOcoQQVKv3aMoCVDedHWgg%40mail.gmail.com.

For more options, visit https://groups.google.com/d/optout.

--
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/CAK4PXd3L43D5mJ-pHQho9AoOF59B7%3DZQtSgZXwhQ5KJSCj_Wow%40mail.gmail.com.

Kong, Fande

unread,
May 22, 2018, 11:02:03 AM5/22/18
to moose...@googlegroups.com
It is still an active research area.  ASM may work for a few hundreds of cores, but it could be problem-dependent. If its performance decreases, you can strength the subdomain solvers by increasing the fill-in level. It may be also helpful to increase the overlap when you see the algorithm does not perform well.


Fande,

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

--
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.

--
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.
To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/CAK4PXd3L43D5mJ-pHQho9AoOF59B7%3DZQtSgZXwhQ5KJSCj_Wow%40mail.gmail.com.
For more options, visit https://groups.google.com/d/optout.

--
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/CAHfiw0e_uxEHKYMTsP6uFDVevxcUCEhw2DgU2G9tevDzHkn7VQ%40mail.gmail.com.

Shohei Ogawa

unread,
May 22, 2018, 11:54:29 AM5/22/18
to moose...@googlegroups.com
Fortunately, I am expecting my simulation will run on hundreds cores, so I can try ASM first. I am not sure about the subdomain solvers and fill-in. Is it something I should always explicitly specify? Is there any documents for Moose?

Shohei

--
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.

--
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.

--
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.

--
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.

Kong, Fande

unread,
May 22, 2018, 12:21:26 PM5/22/18
to moose...@googlegroups.com
Fill-in level:

petsc_options_iname = "-sub_pc_factor_levels"
petsc_options_value  = " 0 "

The default value is 0, and you can increase this number (1, 2, 3, 4..) to have a better convergence.  In my experiences, this number is usually between 0 and 5.  Larger number implies more memory usage.  So don't go to a crazy number.

ASM is a domain decomposition method. If you want to learn more details of ASM, the following book is a good resource.

 Barry Smith, Petter Bjorstad, and William Gropp, Domain Decompositions: Parallel Multilevel Methods for Elliptic Partial Differential Equations, Cambridge University Press.


I am going to add a basic introduction to ASM in MOOSE soon.

Fande,

Shohei

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

--
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.

--
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.

--
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.

--
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.

--
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/CAHfiw0eFZPNVb3fT%2Bgb89DVzMAo9dd2Oknyr41SvuYKTfj8v4Q%40mail.gmail.com.

Derek Gaston

unread,
May 22, 2018, 1:36:18 PM5/22/18
to moose...@googlegroups.com
There are Hypre options for helping with less-elliptic equations.  I definitely recommend taking a spin through our Hypre documentation:


And then the actual Hypre documentation... the User's and Reference manuals are here:


But - what Fande is says is true... in general MG methods are not the best for strongly convective equations.

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.

--
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.

Shohei Ogawa

unread,
May 22, 2018, 4:19:21 PM5/22/18
to moose-users
Fande,
I will try the options. The documentation about ASM for Moose would be really great. Also should I specify sub preconditioner when I use ASM too?

Derek,
I have been checking those hypre documents last few days. Do you have any idea which options have more chance for convective problems?

Thank you,
Shohei 

Kong, Fande

unread,
May 22, 2018, 4:29:41 PM5/22/18
to moose...@googlegroups.com
On Tue, May 22, 2018 at 2:19 PM, Shohei Ogawa <ogawa...@gmail.com> wrote:
Fande,
I will try the options. The documentation about ASM for Moose would be really great. Also should I specify sub preconditioner when I use ASM too?

If ASM with the default setting is already working for you. You do not need to do anything. 


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

Shohei Ogawa

unread,
May 22, 2018, 4:43:48 PM5/22/18
to moose-users
ASM worked for a simple 2D model even though it took so many (hundreds) iterations. When I move to my real problem in 3D, the convergence is not available yet, so I was wondering if I can optimize the settings for ASM.

Thank you,
Shohei

Kong, Fande

unread,
May 22, 2018, 5:12:03 PM5/22/18
to moose...@googlegroups.com
Then go ahead to tune the parameters as I stated in the previous email.  Also pay attention to the memory usage.

Fande

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/bac9b567-f4bc-4dad-9a73-15a2f134c454%40googlegroups.com.

Shohei Ogawa

unread,
May 22, 2018, 6:47:53 PM5/22/18
to moose-users
Thank you! I will try some options.

I also noticed that negative solution value during solve can be a cause of bad convergence. Is there any way to constraint variable values to a certain range, e.g. to be positive? Maybe the way to prevent PJFNK from going too far to be negative for some variables? The physics of my problem includes chemical reactions. There are several variables for species concentration. In my residual and jacobian expressions, 0 is returned if those concentration values are less than 1e-16. Also, they include a natural log of concentration values, I have to do something with negative concentration to avoid error. 

Negative concentration is not physically possible and should not be achieved in numerical simulations. As I stated above, my residual and jacobian expression are force to return 0 and may not be accurate when the concentration is negative. I found that after the first non-linear step, the concentration of one of the species at some part of the domain gets negative. So this may be the cause of bad convergence after the first iteration. After so many non-linear iterations and convergence, the concentration value seems to be in the right range.

Shohei

Sebastian Schunert

unread,
May 22, 2018, 7:04:06 PM5/22/18
to moose...@googlegroups.com
We have similar issues in neutronics that the solution of certain types of problems can become negative [eigenvalue problem solution that converge to a non-fundamental mode]. We do not want these solutions
although they satisfy the equations and can therefore be hit by the eigenvalue solver. Long story short I talked to PETSc people, actually Dmitry Karpeyev back then, about options for constraining the variables to be positive and for reasons
I cannot exactly recall anymore he discouraged this approach quite strongly. Aside from that I doubt that there are the hooks in MOOSE and PETSc to reject certain variable states during the JFNK solve or constrain the attainable solution space, but
I don't know enough to say this with certainty. 

Have you considered to use dampers instead to initially constrain the updates in your Newton solve to be smaller?

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/bfea3c03-097d-46ea-9194-7ad06915b44e%40googlegroups.com.

Shohei Ogawa

unread,
May 22, 2018, 7:55:56 PM5/22/18
to moose-users
Using dampers may be solving the issue. I found a damper which can constrain the update of variable values to a certain range although I don't see many people talking about this feature on the mailing list.



Shohei

Shohei Ogawa

unread,
May 22, 2018, 8:35:25 PM5/22/18
to moose...@googlegroups.com
So far, a dumper to constrain the concentration value above zero didn't work. The concentration values are always positive but didn't change much. So non-linear residual didn't decrease at all. Without dampers, the concentration value got negative at the first non-linear iteration, and gradually increased to be positive. Meanwhile, the non-linear residual decreases. This looks strange but somehow worked.

Shohei

Shohei Ogawa

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/7ba104ec-c166-4e85-81e6-85f3a0f74800%40googlegroups.com.

Sebastian Schunert

unread,
May 22, 2018, 9:52:54 PM5/22/18
to moose...@googlegroups.com
I thought of a less restrictive damper that initially reduces the newton update step size and then limits to not constrain later in the solve. 

This damper does actually exactly what I thought didn't exist - for nodal variables anyway. Good to know but not sure how well it could ever work. 


Shohei Ogawa

To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/7ba104ec-c166-4e85-81e6-85f3a0f74800%40googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

--
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/CAHfiw0ctc4i%3DA32HMtNdGupTPiDao6bZ4mnAvvnScQ_Lmy8HUg%40mail.gmail.com.

Shohei Ogawa

unread,
May 23, 2018, 3:21:03 PM5/23/18
to moose-users
I set a constant damper instead of the bounding value damper. The constant value at 0.5 and 0.2 resulted in the negative value at the first non-linear iteration. 0.01 gave a positive value all over the domain, but it seems that the residual decrease is very slow.

I found that PETSc has a solver setting for inequalities.

The PETSc 3.9 document (page 145) says

5.7 Variational Inequalities

SNES can also solve variational inequalities with box constraints. That are nonlinear algebraic systems with additional inequality constraints on some or all of the variables: Lui  <=ui  <=Hui. Some or all of the lower bounds may be negative infinity (indicated to PETSc with SNES_VI_NINF) and some or all of the upper bounds may be infinity (indicated by SNES_VI_INF). The command 

SNESVISetVariableBounds(SNES,Vec Lu,Vec Hu);
is used to indicate that one is solving a variational inequality. The option -snes_vi_monitor turns on extra monitoring of the active set associated with the bounds and -snes_vi_type allows selecting from several VI solvers, the default is preferred. 

I am not sure if we can use this with Moose. In other words, I am not sure if it is possible to call this function via command line options.

 This might be related to what Sebastian said in the previous post.
Long story short I talked to PETSc people, actually Dmitry Karpeyev back then, about options for constraining the variables to be positive and for reasons
I cannot exactly recall anymore he discouraged this approach quite strongly. Aside from that I doubt that there are the hooks in MOOSE and PETSc to reject certain variable states during the JFNK solve or constrain the attainable solution space, but
I don't know enough to say this with certainty.  

Another approach may be modifying residual or Jacobian for the negative concentration value. Currently, in my implementation, those are 0 for the negative concentration value, so once the temporary solution falls into negative, the residual may be small. Also, even though the concentration gets positive after many non-linear iterations, the solution may be very different from what expected.

Thank you,
Shohei

Kong, Fande

unread,
May 23, 2018, 4:04:07 PM5/23/18
to moose...@googlegroups.com
On Wed, May 23, 2018 at 1:21 PM, Shohei Ogawa <ogawa...@gmail.com> wrote:
I set a constant damper instead of the bounding value damper. The constant value at 0.5 and 0.2 resulted in the negative value at the first non-linear iteration. 0.01 gave a positive value all over the domain, but it seems that the residual decrease is very slow.

I found that PETSc has a solver setting for inequalities.

The PETSc 3.9 document (page 145) says

5.7 Variational Inequalities

SNES can also solve variational inequalities with box constraints. That are nonlinear algebraic systems with additional inequality constraints on some or all of the variables: Lui  <=ui  <=Hui. Some or all of the lower bounds may be negative infinity (indicated to PETSc with SNES_VI_NINF) and some or all of the upper bounds may be infinity (indicated by SNES_VI_INF). The command 

SNESVISetVariableBounds(SNES,Vec Lu,Vec Hu);
is used to indicate that one is solving a variational inequality. The option -snes_vi_monitor turns on extra monitoring of the active set associated with the bounds and -snes_vi_type allows selecting from several VI solvers, the default is preferred. 

I am not sure if we can use this with Moose. In other words, I am not sure if it is possible to call this function via command line options.

We actually have this interface in MOOSE. But we need some resources to make this work properly.

However, I personally do NOT encourage you to use this algorithm for your problems. The constrained nonlinear problem is much harder to solve than a regular nonlinear problem.  SNESVI is originally designed for complementary problems such as mechanics contact with a Lagrange multiplier.

In your original post, I did not see any constraint for the concentration in the equations.  If possible, I would say just make the residual and Jacobian mathematically right even though it may not make sense physically. Newton will eventually give you the right solution after it is converged.  It is not easy to force Newton produce the physical solution during the solving process.  Of course, it should produce the right solution at the end of solving.


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

Shohei Ogawa

unread,
May 23, 2018, 5:42:05 PM5/23/18
to moose-users
The constraint comes to my reaction term used in BC. I didn't include this in the previous post. 

As in the slides attached, this BC models the reaction rate of multiple species depending on the local concentration of multiple species and the electric potential of electrolyte and electrode. The reaction rate per unit area j is

where c_O2 and c_H+ are the concentration of oxygen and proton. \eta is the overpotential of the reaction defined as:

\Psi_electrolyte is the electric potential appearing in the Poisson-Nernst-Planck equation.


Because of the log_e function, I am not sure what value should I return for the residual and Jacobian when the concentration is negative. Other than that I have correctly implemented the residual and Jacobian.

Shohei

pnp_equations.pdf

Shohei Ogawa

unread,
May 23, 2018, 5:42:17 PM5/23/18
to moose-users
The constraint comes to my reaction term used in BC. I didn't include this in the previous post. 

As in the slides attached, this BC models the reaction rate of multiple species depending on the local concentration of multiple species and the electric potential of electrolyte and electrode. The reaction rate per unit area j is

where c_O2 and c_H+ are the concentration of oxygen and proton. \eta is the overpotential of the reaction defined as:

\Psi_electrolyte is the electric potential appearing in the Poisson-Nernst-Planck equation.


Because of the log_e function, I am not sure what value should I return for the residual and Jacobian when the concentration is negative. Other than that I have correctly implemented the residual and Jacobian.

Shohei


pnp_equations.pdf

Sebastian Schunert

unread,
May 23, 2018, 6:14:01 PM5/23/18
to moose...@googlegroups.com
> I set a constant damper instead of the bounding value damper. The constant value at 0.5 and 0.2 resulted in the negative value at the first non-linear iteration. 0.01 gave a positive value all over the domain, but it seems that the residual decrease is very slow.

That is expected since the damper only allows you 1% of the step you'd take without the damper. If nothing else works for you, you can try implementing a custom damper that starts with a small damping parameter and then relaxes it to 1 as the danger of negative solutions becomes smaller.

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/d5bbfb1c-2830-42a9-856a-bdc119d210bc%40googlegroups.com.

Derek Gaston

unread,
May 23, 2018, 6:31:31 PM5/23/18
to moose...@googlegroups.com
Dampers are the right idea here.  This isn't quite the same as the situation Sebastian is in... your _converged_ solution does not have negative values so you just need to suppress them during the nonlnear solve.

Were you using BoundingValueNodalDamper ?  If so... I'm surprised it didn't work well for you because it allows the largest step possible that stays within the bounds.

Are you starting with concentrations that are _zero_?  If so... that's why it isn't moving... because there's no way to allow the solution to move from zero if it's trying to go down.

As Sebastian says: you may need a custom damper... but I really think BoundingValueNodalDamper should work...

Derek

To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/d5bbfb1c-2830-42a9-856a-bdc119d210bc%40googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

--
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/CACOt%3DcUKZGjXhn8VyLBOikW%2Bz_x0UeCOwieq9huz3hv7BDKRRQ%40mail.gmail.com.

Shohei Ogawa

unread,
May 23, 2018, 6:47:07 PM5/23/18
to moose...@googlegroups.com
I started with the concentration value above zero. Since I have the same value for dirichlet BC, the concentration in the domain should be between the initial (boundary) value and zero around the reacting surface. The reaction is consuming those species. This is achieved at the end of non-linear iterations even though the solution goes to negative first at the beginning.

 I have to check how the solution behaves during non linear iterations with BoundingValueNodalDamper and why the residual doesn't decrease.

Shohei

Wang, Yaqi

unread,
May 23, 2018, 7:09:16 PM5/23/18
to moose-users
Not sure if you are aware of this: http://www.mcs.anl.gov/petsc/documentation/faq.html -snes_linesearch_monitor  could be helpful to you.

Shohei

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

--
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.

--
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.

--
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.
To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/CAHfiw0eeHsyCd%3D7MSL4Cs7VMk%2Bj5%2BqOXp6FFG7R4Ch3x5H5SHg%40mail.gmail.com.

Shohei Ogawa

unread,
May 23, 2018, 8:40:26 PM5/23/18
to moose...@googlegroups.com
It seems that BoundingValueNodalDamper doesn't work because after a few non-linear iterations the damping factor gets 0. So the solution doesn't move at all.

Does this block seem reasonable? I am not sure what value I should put for min_damping.

[c_O2_positive_damper]

type = BoundingValueNodalDamper

min_value = 1e-16

variable = c_O2

max_value = 1e10

min_damping = -1e4

[]


Also, should I specify specific line search methods?


With line_search = 'bt' for example, the concentration actually got negative even with the dumper.


Thank you,

Shohei



Shohei Ogawa

To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/CABk6%3D%2BbEpAUohS8V6HB8FgmaiNMKDLVHhx7067aXjxhnUwMsdw%40mail.gmail.com.

Shohei Ogawa

unread,
Jun 13, 2018, 5:53:20 PM6/13/18
to moose...@googlegroups.com
I just want to report back here. I somehow found the cause of the issue, where the concentration solution goes negative in the first non-linear iteration and converged result seems to be unphysical. I was using a linear Lagrange approximation for the species concentration variables and a second order Lagrange approximation for the electric potential variable where its gradient drives the movement of charged ion. This caused the issue. After I changed the order for the electric potential variable, from the second to first, the simulation works well and the convergence is much better.

Thank you,
Shohei

Shohei Ogawa
Reply all
Reply to author
Forward
0 new messages