Re: [Libmesh-users] Jacobian for Nonlinear Solve

 Re: [Libmesh-users] Jacobian for Nonlinear Solve From: John Peterson - 2008-07-23 14:10:31 ```On Wed, Jul 23, 2008 at 8:52 AM, Derek Gaston wrote: > > Sorry for the long winded reply... the short of it is: just use a global > variable ;-) Ouch. Yeah this needs to not be the "official" way of doing this. Unless there is a technical problem with changing the signature of the jacobian and residual functions (?) I would second the notion of having it take (probably a pointer to, so it could be NULL by default) an EquationSystems object. -- John ```

Thread view

 [Libmesh-users] Jacobian for Nonlinear Solve From: Andrea Hawkins - 2008-07-22 23:30:03 ```Hi- I am trying to solve a nonlinear system using as my system class TransientNonlinearImplicitSystem with the PetscNonlinearSolver. I realize the PetscNonlinearSolver class has a function pointer for jacobian, but I am wondering how when the function jacobian is actually defined to interface with the TransientNonlinearImplicitSystem. i.e. how do I access the equation system. Do you have a simple example of something similar being done I could look at? Thanks, Andrea ```
 Re: [Libmesh-users] Jacobian for Nonlinear Solve From: John Peterson - 2008-07-23 00:32:05 ```Hi, On Tue, Jul 22, 2008 at 6:30 PM, Andrea Hawkins wrote: > Hi- > > I am trying to solve a nonlinear system using as my system class > TransientNonlinearImplicitSystem with the PetscNonlinearSolver. I > realize the PetscNonlinearSolver class has a function pointer for > jacobian, but I am wondering how when the function jacobian is > actually defined to interface with the > TransientNonlinearImplicitSystem. i.e. how do I access the equation > system. A System can return a reference to the EquationSystem of which it is a part. system.get_equation_systems(); should do it (see include/solvers/system.h). Does that answer your question? > Do you have a simple example of something similar being done I could look at? I believe Derek's writing nonlinear solvers right now, so he will probably have more to say shortly... -- John ```
 Re: [Libmesh-users] Jacobian for Nonlinear Solve From: Andrea Hawkins - 2008-07-23 02:19:19 ```This must be a syntax error that I'm not seeing, but when I try the command EquationSystems& es = system.get_equation_systems(); in the jacobian function, I get the error error: request for member get_equation_systems in system, which is of non-class type. Thanks, Andrea _ It seems like this function might not have access to it in this manner, since only the function pointer is in the class nonlinear_solver and not the function itself. On Tue, Jul 22, 2008 at 7:32 PM, John Peterson wrote: > Hi, > > On Tue, Jul 22, 2008 at 6:30 PM, Andrea Hawkins wrote: >> Hi- >> >> I am trying to solve a nonlinear system using as my system class >> TransientNonlinearImplicitSystem with the PetscNonlinearSolver. I >> realize the PetscNonlinearSolver class has a function pointer for >> jacobian, but I am wondering how when the function jacobian is >> actually defined to interface with the >> TransientNonlinearImplicitSystem. i.e. how do I access the equation >> system. > > > A System can return a reference to the EquationSystem of which it is a part. > > system.get_equation_systems(); > > should do it (see include/solvers/system.h). Does that answer your question? > > >> Do you have a simple example of something similar being done I could look at? > > I believe Derek's writing nonlinear solvers right now, so he will > probably have more to say shortly... > > -- > John > ```
 Re: [Libmesh-users] Jacobian for Nonlinear Solve From: John Peterson - 2008-07-23 11:51:50 ```On Tue, Jul 22, 2008 at 9:19 PM, Andrea Hawkins wrote: > This must be a syntax error that I'm not seeing, but when I try the command > > EquationSystems& es = system.get_equation_systems(); > > in the jacobian function, I get the error > > error: request for member get_equation_systems in system, which is of > non-class type. >From within a NonlinearSolver or PetscNonlinearSolver, you will need to do: EquationSystems& es = system().get_equation_systems(); since NonlinearSolver::system() is a member function returning a reference to the present System. My previous email was a little confusing, I was just trying to convey the idea that get_equation_systems() was a System:: member. -- John ```
 Re: [Libmesh-users] Jacobian for Nonlinear Solve From: John Peterson - 2008-07-23 12:18:57 ```On Tue, Jul 22, 2008 at 9:19 PM, Andrea Hawkins wrote: > > It seems like this function might not have access to it in this > manner, since only the function pointer is in the class > nonlinear_solver and not the function itself. I think I finally understand what the problem is... the "jacobian" is actually a function pointer that you attach, and therefore it knows nothing about the contents of the NonlinearSolver class, only the Matrix and Vector that are passed in. Hm...I'm not immediately sure how to get around this issue, but I'll keep looking. -- John ```
 Re: [Libmesh-users] Jacobian for Nonlinear Solve From: Derek Gaston - 2008-07-23 13:53:02 Attachments: Message as HTML ```Check out src/numerics/petsc_nonlinear_solver.C... around line 112 there is a __petsc_snes_jacobian function. This is the _actual_ function that libMesh passes into Petsc SNES for you. The genius is that you can pass one pointer (which comes through as "ctx" in __petsc_snes_jacobian() ) to Petsc that Petsc will just pass back to you all the time. libMesh (Ben?) takes advantage of this by passing a pointer to the solver. That means that "ctx" is always a pointer to the solver... but to use it we use a static_cast to get it back to a solver pointer. From there, you can get to the system (as it does there in __petsc_snes_jacobian. Now... that still might not answer your question because the jacobian assembly function that you attach to the nonlinear-implicit-system has a signature of: (const NumericVector&, SparseMatrix&) which obviously doesn't pass the system. So what to do? In Ben's example code that uses NonlinearImplicitSystem (I've attached it just in case you missed it a few months ago) he just does this at the top of his main .C file: EquationSystems *_equation_system = NULL; essentially just creating a global variable that holds a pointer to the equation system. He sets that in his main... and then his compute_jacobian() function just uses the global variable.... Personally, I don't use that approach... but mine is even more complicated (I have objects that provide residuals / jacobians and they store references to their own EquationSystems object...). I don't see why we couldn't modify the signature of the function pointer that you attach to NonlinearImplicitSystem to also pass an EquationSystems object. It's kind of hard to assemble anything without it.... and we're already getting pointers to it in __petsc_snes_jacobian Anyone have any objections? Sorry for the long winded reply... the short of it is: just use a global variable ;-) Derek PS - Ben... we really need to make this an official example. Do you want me to do it? ```
 Re: [Libmesh-users] Jacobian for Nonlinear Solve From: John Peterson - 2008-07-23 14:10:31 ```On Wed, Jul 23, 2008 at 8:52 AM, Derek Gaston wrote: > > Sorry for the long winded reply... the short of it is: just use a global > variable ;-) Ouch. Yeah this needs to not be the "official" way of doing this. Unless there is a technical problem with changing the signature of the jacobian and residual functions (?) I would second the notion of having it take (probably a pointer to, so it could be NULL by default) an EquationSystems object. -- John ```
 Re: [Libmesh-users] Jacobian for Nonlinear Solve From: Roy Stogner - 2008-07-23 14:24:06 ``` On Wed, 23 Jul 2008, John Peterson wrote: > On Wed, Jul 23, 2008 at 8:52 AM, Derek Gaston wrote: >> >> Sorry for the long winded reply... the short of it is: just use a global >> variable ;-) > > Ouch. Yeah this needs to not be the "official" way of doing this. > Unless there is a technical problem with changing the signature of the > jacobian and residual functions (?) I would second the notion of > having it take (probably a pointer to, so it could be NULL by default) > an EquationSystems object. In that case you might still have to have some global to "remember" which system in the EquationSystems you were solving, if you had two loosely coupled systems and wanted to use PETSc on both. Better to pass a pointer to a System object, and get the EquationSystems object from that when necessary. But anyway, I agree with John's main point: we ought to be passing more information in to *residual/*jacobian/*matvec. Using the existing interface is awkward enough that I suspect even those people whose code will be broken by an API change will thank us. --- Roy ```
 Re: [Libmesh-users] Jacobian for Nonlinear Solve From: Andrea Hawkins - 2008-07-23 14:50:32 ```All right, so for now I will just use a global variable. =) I am wondering though if it will work to just change the signature of the jacobian function because when Petsc calls it, it needs the current signature. (at least, that's how i understand it...) Thanks for the help! Andrea On Wed, Jul 23, 2008 at 9:23 AM, Roy Stogner wrote: > > > On Wed, 23 Jul 2008, John Peterson wrote: > >> On Wed, Jul 23, 2008 at 8:52 AM, Derek Gaston wrote: >>> >>> Sorry for the long winded reply... the short of it is: just use a global >>> variable ;-) >> >> Ouch. Yeah this needs to not be the "official" way of doing this. >> Unless there is a technical problem with changing the signature of the >> jacobian and residual functions (?) I would second the notion of >> having it take (probably a pointer to, so it could be NULL by default) >> an EquationSystems object. > > In that case you might still have to have some global to "remember" > which system in the EquationSystems you were solving, if you had two > loosely coupled systems and wanted to use PETSc on both. Better to > pass a pointer to a System object, and get the EquationSystems object > from that when necessary. > > But anyway, I agree with John's main point: we ought to be passing > more information in to *residual/*jacobian/*matvec. Using the > existing interface is awkward enough that I suspect even those people > whose code will be broken by an API change will thank us. > --- > Roy > > ------------------------------------------------------------------------- > This SF.Net email is sponsored by the Moblin Your Move Developer's challenge > Build the coolest Linux based applications with Moblin SDK & win great prizes > Grand prize is a trip for two to an Open Source event anywhere in the world > http://moblin-contest.org/redirect.php?banner_id=100&url=/ > _______________________________________________ > Libmesh-users mailing list > Libmesh-users@... > https://lists.sourceforge.net/lists/listinfo/libmesh-users > ```
 Re: [Libmesh-users] Jacobian for Nonlinear Solve From: Roy Stogner - 2008-07-23 15:10:05 ```On Wed, 23 Jul 2008, Andrea Hawkins wrote: > I am wondering though if it will work to just change the signature of > the jacobian function because when Petsc calls it, it needs the > current signature. (at least, that's how i understand it...) Petsc needs the current signature of __libmesh_petsc_snes_residual and the other global functions we define in petsc_nonlinear_solver.C, but we're free to change around the signatures of the user-specified functions they call. We've even got a little leeway with the __libmesh_petsc_ functions, since we could decide to use the catch-all "void *ctx" in a different way. The trouble with changing around PetscNonlinearSolver::jacobian and the like isn't compatibility with PETSc code, it's compatibility with user code. I warned you before that we'd eventually have to change the DiffSystem assembly APIs to work better with threading; now it turns out we're going to have to change the non-DiffSystem nonlinear solver assembly APIs too? I think both changes are worth the annoyance, but (along with Mesh/MeshBase, prefixes in asserts and Makefile constants, etc.) it feels like we're accumulating too many little API changes for our next release. People using 0.6.2 may be upset at all the little thing they'll have to fiddle with to move to 0.6.9. --- Roy ```