Hi everyone, I am attempting to solve the electrolyte domain of li-ion battery model. The concentration problem is somewhat like a nonlinear diffusion/heat equation but it gets messy because it has a div(grad(ln(c))) term. So I would be looking at something like this: -div(k1*grad(c)) + div(k2*grad(ln(c))) = 0 Which we can rewrite as: -div((k1 - k2/c)*grad(c)) = 0 If k2 is relatively small enough, the newton solver should converge with little problem. The problem I am having is with the linear solve. If I discretize the above using CG1, no solver works so long as k2 >> 0. With DG1, I can only use direct solvers like LU to solve the resulting linear system of equations. I can't seem to find a good iterative solver/pc to solve this - HYPRE gives me a DIVERGED_NANORINF, and jacobi/bjacobi do not converge at all. I've attached the code for you guys to play around with. I don't know what is a good PC for grad(ln(c)), my guess is that I'm going to need to provide an assembled preconditioner which approximates the inverse of the k2/c*grad(c) term but it doesn't seem clear to me how I can provide an explicit PC to the nonlinear solver without resorting to hacking at petsc4py. Any thoughts appreciated Justin
On 7 Nov 2018, at 08:30, Justin Chang <jychang48@gmail.com> wrote:
Hi everyone,
I am attempting to solve the electrolyte domain of li-ion battery model. The concentration problem is somewhat like a nonlinear diffusion/heat equation but it gets messy because it has a div(grad(ln(c))) term. So I would be looking at something like this:
-div(k1*grad(c)) + div(k2*grad(ln(c))) = 0
Which we can rewrite as:
-div((k1 - k2/c)*grad(c)) = 0
If k2 is relatively small enough, the newton solver should converge with little problem. The problem I am having is with the linear solve. If I discretize the above using CG1, no solver works so long as k2 >> 0. With DG1, I can only use direct solvers like LU to solve the resulting linear system of equations. I can't seem to find a good iterative solver/pc to solve this - HYPRE gives me a DIVERGED_NANORINF, and jacobi/bjacobi do not converge at all.
Not sure about preconditioning immediately. But how are you ensuring that you never have c=0 in your solves?
I've attached the code for you guys to play around with. I don't know what is a good PC for grad(ln(c)), my guess is that I'm going to need to provide an assembled preconditioner which approximates the inverse of the k2/c*grad(c) term but it doesn't seem clear to me how I can provide an explicit PC to the nonlinear solver without resorting to hacking at petsc4py.
You can use an ExplicitSchurPC (see eg usage here https://firedrakeproject.org/demos/geometric_multigrid.py.html) along with pc typ python. No need to write petsc4py code directly. Lawrence
Justin
<nonlinear_poisson_pc.py> _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
On Wed, Nov 7, 2018 at 4:02 AM Lawrence Mitchell <wence@gmx.li> wrote:
On 7 Nov 2018, at 08:30, Justin Chang <jychang48@gmail.com> wrote:
Hi everyone,
I am attempting to solve the electrolyte domain of li-ion battery model. The concentration problem is somewhat like a nonlinear diffusion/heat equation but it gets messy because it has a div(grad(ln(c))) term. So I would be looking at something like this:
-div(k1*grad(c)) + div(k2*grad(ln(c))) = 0
Which we can rewrite as:
-div((k1 - k2/c)*grad(c)) = 0
If k2 is relatively small enough, the newton solver should converge with little problem. The problem I am having is with the linear solve. If I discretize the above using CG1, no solver works so long as k2 >> 0. With DG1, I can only use direct solvers like LU to solve the resulting linear system of equations. I can't seem to find a good iterative solver/pc to solve this - HYPRE gives me a DIVERGED_NANORINF, and jacobi/bjacobi do not converge at all.
Not sure about preconditioning immediately. But how are you ensuring that you never have c=0 in your solves?
I've attached the code for you guys to play around with. I don't know what is a good PC for grad(ln(c)), my guess is that I'm going to need to provide an assembled preconditioner which approximates the inverse of the k2/c*grad(c) term but it doesn't seem clear to me how I can provide an explicit PC to the nonlinear solver without resorting to hacking at petsc4py.
You can use an ExplicitSchurPC (see eg usage here https://firedrakeproject.org/demos/geometric_multigrid.py.html) along with pc typ python. No need to write petsc4py code directly.
I must be reading the code wrong because I do not see the coefficient. For the system in the file, Hypre should work fine. More broadly, does the cparse grid solve converge? If so, can you show the smoother iterates for a 2-level problem? Thanks, Matt
Lawrence
Justin
<nonlinear_poisson_pc.py>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
On 7 Nov 2018, at 09:56, Matthew Knepley <knepley@gmail.com> wrote:
On Wed, Nov 7, 2018 at 4:02 AM Lawrence Mitchell <wence@gmx.li> wrote:
On 7 Nov 2018, at 08:30, Justin Chang <jychang48@gmail.com> wrote:
Hi everyone,
I am attempting to solve the electrolyte domain of li-ion battery model. The concentration problem is somewhat like a nonlinear diffusion/heat equation but it gets messy because it has a div(grad(ln(c))) term. So I would be looking at something like this:
-div(k1*grad(c)) + div(k2*grad(ln(c))) = 0
Which we can rewrite as:
-div((k1 - k2/c)*grad(c)) = 0
If k2 is relatively small enough, the newton solver should converge with little problem. The problem I am having is with the linear solve. If I discretize the above using CG1, no solver works so long as k2 >> 0. With DG1, I can only use direct solvers like LU to solve the resulting linear system of equations. I can't seem to find a good iterative solver/pc to solve this - HYPRE gives me a DIVERGED_NANORINF, and jacobi/bjacobi do not converge at all.
Not sure about preconditioning immediately. But how are you ensuring that you never have c=0 in your solves?
I've attached the code for you guys to play around with. I don't know what is a good PC for grad(ln(c)), my guess is that I'm going to need to provide an assembled preconditioner which approximates the inverse of the k2/c*grad(c) term but it doesn't seem clear to me how I can provide an explicit PC to the nonlinear solver without resorting to hacking at petsc4py.
You can use an ExplicitSchurPC (see eg usage here https://firedrakeproject.org/demos/geometric_multigrid.py.html) along with pc typ python. No need to write petsc4py code directly.
I must be reading the code wrong because I do not see the coefficient. For the system in the file, Hypre should work fine. More broadly, does the cparse grid solve converge? If so, can you show the smoother iterates for a 2-level problem?
If we ignore the IP scheme and just use the C0 scheme Justin posted: Here's the residual, with u_ the state to be solved for: F = dot(grad(v), D(u_)*grad(u_))*dx - v*f*dx - g*v*ds(3) - g*v*ds(4) and D is defined as: def D(c): return Constant(100) - Constant(1.0)/c Which gives you the coefficient variation. I worry with this, given that the linearisation about u* looks like: grad(v) : (grad(u) (k1 - k2/u*) + grad(u*) k2 u / u*^2) So if u* is zero, this blows up, no? Cheers, Lawrence
On Wed, Nov 7, 2018 at 5:24 AM Lawrence Mitchell <wence@gmx.li> wrote:
On 7 Nov 2018, at 09:56, Matthew Knepley <knepley@gmail.com> wrote:
On Wed, Nov 7, 2018 at 4:02 AM Lawrence Mitchell <wence@gmx.li> wrote:
On 7 Nov 2018, at 08:30, Justin Chang <jychang48@gmail.com> wrote:
Hi everyone,
I am attempting to solve the electrolyte domain of li-ion battery model. The concentration problem is somewhat like a nonlinear diffusion/heat equation but it gets messy because it has a div(grad(ln(c))) term. So I would be looking at something like this:
-div(k1*grad(c)) + div(k2*grad(ln(c))) = 0
Which we can rewrite as:
-div((k1 - k2/c)*grad(c)) = 0
If k2 is relatively small enough, the newton solver should converge with little problem. The problem I am having is with the linear solve. If I discretize the above using CG1, no solver works so long as k2 >> 0. With DG1, I can only use direct solvers like LU to solve the resulting linear system of equations. I can't seem to find a good iterative solver/pc to solve this - HYPRE gives me a DIVERGED_NANORINF, and jacobi/bjacobi do not converge at all.
Not sure about preconditioning immediately. But how are you ensuring that you never have c=0 in your solves?
I've attached the code for you guys to play around with. I don't know what is a good PC for grad(ln(c)), my guess is that I'm going to need to provide an assembled preconditioner which approximates the inverse of the k2/c*grad(c) term but it doesn't seem clear to me how I can provide an explicit PC to the nonlinear solver without resorting to hacking at petsc4py.
You can use an ExplicitSchurPC (see eg usage here https://firedrakeproject.org/demos/geometric_multigrid.py.html) along with pc typ python. No need to write petsc4py code directly.
I must be reading the code wrong because I do not see the coefficient. For the system in the file, Hypre should work fine. More broadly, does the cparse grid solve converge? If so, can you show the smoother iterates for a 2-level problem?
If we ignore the IP scheme and just use the C0 scheme Justin posted:
Here's the residual, with u_ the state to be solved for:
F = dot(grad(v), D(u_)*grad(u_))*dx - v*f*dx - g*v*ds(3) - g*v*ds(4)
and D is defined as:
def D(c): return Constant(100) - Constant(1.0)/c
Which gives you the coefficient variation. I worry with this, given that the linearisation about u* looks like:
grad(v) : (grad(u) (k1 - k2/u*) + grad(u*) k2 u / u*^2)
So if u* is zero, this blows up, no?
I see. Yes, you are right. This is similar to the p-Laplacian, except its u in the denominator instead of grad u. There you usually regularize with (u + \epsilon), which should probably be done here. Sometimes you can get away with solving the coarse grid with svd and hope nothing goes wrong in the smoother. Matt
Cheers,
Lawrence
-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
Lawrence, the short answer is no the concentration in the electrolyte won't reach zero for the problems we've considered because it mainly acts as a "transmitter" of Li-ions between the two electrode interfaces. We're more worried about concentrations going below zero in the anode/cathode both those two domains are simple poisson solves. But yes I am worried that in theory the newton iterates could approach zero which is where I'm thinking VI can help prevent that. I'll have a look into the pc_type python stuff. Matt, can you elaborate a little more on how one would construct a PC for problems like the p-Laplacian? I was also thinking something similar should be done since the k2/c term of my problem is clearly screwing up the out-of-box hypre solve. On Wed, Nov 7, 2018 at 4:23 AM Matthew Knepley <knepley@gmail.com> wrote:
On Wed, Nov 7, 2018 at 5:24 AM Lawrence Mitchell <wence@gmx.li> wrote:
On 7 Nov 2018, at 09:56, Matthew Knepley <knepley@gmail.com> wrote:
On Wed, Nov 7, 2018 at 4:02 AM Lawrence Mitchell <wence@gmx.li> wrote:
On 7 Nov 2018, at 08:30, Justin Chang <jychang48@gmail.com> wrote:
Hi everyone,
I am attempting to solve the electrolyte domain of li-ion battery model. The concentration problem is somewhat like a nonlinear diffusion/heat equation but it gets messy because it has a div(grad(ln(c))) term. So I would be looking at something like this:
-div(k1*grad(c)) + div(k2*grad(ln(c))) = 0
Which we can rewrite as:
-div((k1 - k2/c)*grad(c)) = 0
If k2 is relatively small enough, the newton solver should converge with little problem. The problem I am having is with the linear solve. If I discretize the above using CG1, no solver works so long as k2 >> 0. With DG1, I can only use direct solvers like LU to solve the resulting linear system of equations. I can't seem to find a good iterative solver/pc to solve this - HYPRE gives me a DIVERGED_NANORINF, and jacobi/bjacobi do not converge at all.
Not sure about preconditioning immediately. But how are you ensuring that you never have c=0 in your solves?
I've attached the code for you guys to play around with. I don't know what is a good PC for grad(ln(c)), my guess is that I'm going to need to provide an assembled preconditioner which approximates the inverse of the k2/c*grad(c) term but it doesn't seem clear to me how I can provide an explicit PC to the nonlinear solver without resorting to hacking at petsc4py.
You can use an ExplicitSchurPC (see eg usage here https://firedrakeproject.org/demos/geometric_multigrid.py.html) along with pc typ python. No need to write petsc4py code directly.
I must be reading the code wrong because I do not see the coefficient. For the system in the file, Hypre should work fine. More broadly, does the cparse grid solve converge? If so, can you show the smoother iterates for a 2-level problem?
If we ignore the IP scheme and just use the C0 scheme Justin posted:
Here's the residual, with u_ the state to be solved for:
F = dot(grad(v), D(u_)*grad(u_))*dx - v*f*dx - g*v*ds(3) - g*v*ds(4)
and D is defined as:
def D(c): return Constant(100) - Constant(1.0)/c
Which gives you the coefficient variation. I worry with this, given that the linearisation about u* looks like:
grad(v) : (grad(u) (k1 - k2/u*) + grad(u*) k2 u / u*^2)
So if u* is zero, this blows up, no?
I see. Yes, you are right. This is similar to the p-Laplacian, except its u in the denominator instead of grad u. There you usually regularize with (u + \epsilon), which should probably be done here. Sometimes you can get away with solving the coarse grid with svd and hope nothing goes wrong in the smoother.
Matt
Cheers,
Lawrence
-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener
https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
On Wed, Nov 7, 2018 at 11:38 AM Justin Chang <jychang48@gmail.com> wrote:
Lawrence, the short answer is no the concentration in the electrolyte won't reach zero for the problems we've considered because it mainly acts as a "transmitter" of Li-ions between the two electrode interfaces. We're more worried about concentrations going below zero in the anode/cathode both those two domains are simple poisson solves. But yes I am worried that in theory the newton iterates could approach zero which is where I'm thinking VI can help prevent that.
I'll have a look into the pc_type python stuff.
Matt, can you elaborate a little more on how one would construct a PC for problems like the p-Laplacian? I was also thinking something similar should be done since the k2/c term of my problem is clearly screwing up the out-of-box hypre solve.
I have a p-Laplacian example in SNES ex12. Its just GMG, which I think should work here if c is not too crazy. Matt
On Wed, Nov 7, 2018 at 4:23 AM Matthew Knepley <knepley@gmail.com> wrote:
On Wed, Nov 7, 2018 at 5:24 AM Lawrence Mitchell <wence@gmx.li> wrote:
On 7 Nov 2018, at 09:56, Matthew Knepley <knepley@gmail.com> wrote:
On Wed, Nov 7, 2018 at 4:02 AM Lawrence Mitchell <wence@gmx.li> wrote:
On 7 Nov 2018, at 08:30, Justin Chang <jychang48@gmail.com> wrote:
Hi everyone,
I am attempting to solve the electrolyte domain of li-ion battery model. The concentration problem is somewhat like a nonlinear diffusion/heat equation but it gets messy because it has a div(grad(ln(c))) term. So I would be looking at something like this:
-div(k1*grad(c)) + div(k2*grad(ln(c))) = 0
Which we can rewrite as:
-div((k1 - k2/c)*grad(c)) = 0
If k2 is relatively small enough, the newton solver should converge with little problem. The problem I am having is with the linear solve. If I discretize the above using CG1, no solver works so long as k2 >> 0. With DG1, I can only use direct solvers like LU to solve the resulting linear system of equations. I can't seem to find a good iterative solver/pc to solve this - HYPRE gives me a DIVERGED_NANORINF, and jacobi/bjacobi do not converge at all.
Not sure about preconditioning immediately. But how are you ensuring that you never have c=0 in your solves?
I've attached the code for you guys to play around with. I don't know what is a good PC for grad(ln(c)), my guess is that I'm going to need to provide an assembled preconditioner which approximates the inverse of the k2/c*grad(c) term but it doesn't seem clear to me how I can provide an explicit PC to the nonlinear solver without resorting to hacking at petsc4py.
You can use an ExplicitSchurPC (see eg usage here https://firedrakeproject.org/demos/geometric_multigrid.py.html) along with pc typ python. No need to write petsc4py code directly.
I must be reading the code wrong because I do not see the coefficient. For the system in the file, Hypre should work fine. More broadly, does the cparse grid solve converge? If so, can you show the smoother iterates for a 2-level problem?
If we ignore the IP scheme and just use the C0 scheme Justin posted:
Here's the residual, with u_ the state to be solved for:
F = dot(grad(v), D(u_)*grad(u_))*dx - v*f*dx - g*v*ds(3) - g*v*ds(4)
and D is defined as:
def D(c): return Constant(100) - Constant(1.0)/c
Which gives you the coefficient variation. I worry with this, given that the linearisation about u* looks like:
grad(v) : (grad(u) (k1 - k2/u*) + grad(u*) k2 u / u*^2)
So if u* is zero, this blows up, no?
I see. Yes, you are right. This is similar to the p-Laplacian, except its u in the denominator instead of grad u. There you usually regularize with (u + \epsilon), which should probably be done here. Sometimes you can get away with solving the coarse grid with svd and hope nothing goes wrong in the smoother.
Matt
Cheers,
Lawrence
-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener
https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
Okay so I had a slight error in my code. I forgot to transpose the u0 expression by some value such that it won't be negative anywhere (I thought I did it earlier, but I guess I forgot when I did the FEniCS to Firedrake conversion for the MWE). It fixed my issue with CG1, but the IP formulation still does not work. On Wed, Nov 7, 2018 at 10:25 AM Matthew Knepley <knepley@gmail.com> wrote:
On Wed, Nov 7, 2018 at 11:38 AM Justin Chang <jychang48@gmail.com> wrote:
Lawrence, the short answer is no the concentration in the electrolyte won't reach zero for the problems we've considered because it mainly acts as a "transmitter" of Li-ions between the two electrode interfaces. We're more worried about concentrations going below zero in the anode/cathode both those two domains are simple poisson solves. But yes I am worried that in theory the newton iterates could approach zero which is where I'm thinking VI can help prevent that.
I'll have a look into the pc_type python stuff.
Matt, can you elaborate a little more on how one would construct a PC for problems like the p-Laplacian? I was also thinking something similar should be done since the k2/c term of my problem is clearly screwing up the out-of-box hypre solve.
I have a p-Laplacian example in SNES ex12. Its just GMG, which I think should work here if c is not too crazy.
Matt
On Wed, Nov 7, 2018 at 4:23 AM Matthew Knepley <knepley@gmail.com> wrote:
On Wed, Nov 7, 2018 at 5:24 AM Lawrence Mitchell <wence@gmx.li> wrote:
On 7 Nov 2018, at 09:56, Matthew Knepley <knepley@gmail.com> wrote:
On Wed, Nov 7, 2018 at 4:02 AM Lawrence Mitchell <wence@gmx.li> wrote:
On 7 Nov 2018, at 08:30, Justin Chang <jychang48@gmail.com> wrote:
Hi everyone,
I am attempting to solve the electrolyte domain of li-ion battery model. The concentration problem is somewhat like a nonlinear diffusion/heat equation but it gets messy because it has a div(grad(ln(c))) term. So I would be looking at something like this:
-div(k1*grad(c)) + div(k2*grad(ln(c))) = 0
Which we can rewrite as:
-div((k1 - k2/c)*grad(c)) = 0
If k2 is relatively small enough, the newton solver should converge with little problem. The problem I am having is with the linear solve. If I discretize the above using CG1, no solver works so long as k2 >> 0. With DG1, I can only use direct solvers like LU to solve the resulting linear system of equations. I can't seem to find a good iterative solver/pc to solve this - HYPRE gives me a DIVERGED_NANORINF, and jacobi/bjacobi do not converge at all.
Not sure about preconditioning immediately. But how are you ensuring that you never have c=0 in your solves?
I've attached the code for you guys to play around with. I don't know what is a good PC for grad(ln(c)), my guess is that I'm going to need to provide an assembled preconditioner which approximates the inverse of the k2/c*grad(c) term but it doesn't seem clear to me how I can provide an explicit PC to the nonlinear solver without resorting to hacking at petsc4py.
You can use an ExplicitSchurPC (see eg usage here https://firedrakeproject.org/demos/geometric_multigrid.py.html) along with pc typ python. No need to write petsc4py code directly.
I must be reading the code wrong because I do not see the coefficient. For the system in the file, Hypre should work fine. More broadly, does the cparse grid solve converge? If so, can you show the smoother iterates for a 2-level problem?
If we ignore the IP scheme and just use the C0 scheme Justin posted:
Here's the residual, with u_ the state to be solved for:
F = dot(grad(v), D(u_)*grad(u_))*dx - v*f*dx - g*v*ds(3) - g*v*ds(4)
and D is defined as:
def D(c): return Constant(100) - Constant(1.0)/c
Which gives you the coefficient variation. I worry with this, given that the linearisation about u* looks like:
grad(v) : (grad(u) (k1 - k2/u*) + grad(u*) k2 u / u*^2)
So if u* is zero, this blows up, no?
I see. Yes, you are right. This is similar to the p-Laplacian, except its u in the denominator instead of grad u. There you usually regularize with (u + \epsilon), which should probably be done here. Sometimes you can get away with solving the coarse grid with svd and hope nothing goes wrong in the smoother.
Matt
Cheers,
Lawrence
-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener
https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener
https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
On 07/11/2018 19:52, Justin Chang wrote:
Okay so I had a slight error in my code. I forgot to transpose the u0 expression by some value such that it won't be negative anywhere (I thought I did it earlier, but I guess I forgot when I did the FEniCS to Firedrake conversion for the MWE). It fixed my issue with CG1, but the IP formulation still does not work.
Having a quick look at your code, I think the penalty term in your IP formulation: + alpha/h_avg*dot(jump(v, n), jump(u_, n))*dS \ needs scaling by D_u. Ideally there should be a factor corresponding to the local min/max ratio of D_u as well. Have a look at http://dx.doi.org/10.1016/j.cam.2006.08.029 for reference. In general, I don't think IP works very well for strongly varied diffusion coefficients, the penalty is either too small and things go unstable or too large and the system becomes very stiff Cheers Stephan
participants (4)
-
Justin Chang
-
Lawrence Mitchell
-
Matthew Knepley
-
Stephan Kramer