From: Dmitry Karpeev <karpeev@mc...>  20120509 23:14:45

Ata, system.solve() will call it. I assume that by "changing constantly" you mean "once per timestep"? SNESSolve() initiates the call that eventually dispatches to your bounds calculation routine. SNES calls that once per solve. Dmitry. On Wed, May 9, 2012 at 5:22 PM, Ataollah Mesgarnejad <amesga1@... > wrote: > Dmitry, > > A quick question: When does the bounds routine get called? My bounds are > changing constantly do I need to call the bounds function myself or does it > get called for each system.solve()? > > Thanks, > Ata > > On May 5, 2012, at 12:21 PM, Dmitry Karpeev wrote: > > Okay, that helps. > I'm looking further into it. > Dmitry. > > On Sat, May 5, 2012 at 12:08 PM, Ataollah Mesgarnejad < > amesga1@...> wrote: > >> Dmitry, >> >> I just tried it without bounds and I get the same DIVERGED_BREAKDOWN >> error. >> >> Ata >> >> Ata, >> I'm trying to reproduce the problem that you are observing, so I wanted >> to understand better what you are doing there. >> Do you just skip Jacobian evaluation and return from the ComputeJacobian >> call immediately? >> Also, without usepetscdm, how are you able to run with VIRS (I >> remember you mentioned your own wrapper to SNESVI)? Finally, do you see >> this problem when running with SNES (not SNESVI)? I wonder if the issue is >> in the interaction with the bounds. >> >> Thanks for your patience. >> Dmitry. >> >> On Wed, May 2, 2012 at 7:51 AM, Ataollah Mesgarnejad < >> amesga1@...> wrote: >> >>> >>> On May 1, 2012, at 10:16 PM, Dmitry Karpeev wrote: >>> >>> Yeah, it's a bit of a hack right now  DMCreateMatrix simply increfs >>> and return the matrix from the underlying NonlinearImplicitSystem. The >>> reason is that initially the matrix is preallocated by unassembled, and >>> MatDuplicate >>> doesn't want to duplicate an unassembled matrix. I can fix that, but >>> I'm a bit puzzled by the fact that KSP messes >>> with the matrix. Can you run with snes_view and send me the output? >>> >>> Dmitry, >>> >>> Here it is: >>> >>>  with the flag  >>> >>> Linear solve converged due to CONVERGED_RTOL iterations 44 >>> Calculating surface energy... done >>> Calculating fracture energy... done >>> elastic energy:13.74054,fracture energy:0.00000,total:13.74054 >>> applying variable bounds ... done >>> NL step 0, residual_2 = 2.41046e01 >>> SNES Jacobian assembly... done >>> Linear solve converged due to CONVERGED_RTOL iterations 3 >>> NL step 1, residual_2 = 5.10144e07 >>> Linear solve did not converge due to DIVERGED_BREAKDOWN iterations 1 >>> SNES Object: 1 MPI processes >>> type: virs >>> maximum iterations=50, maximum function evaluations=10000 >>> tolerances: relative=1e08, absolute=1e35, solution=1e08 >>> total number of linear solver iterations=3 >>> total number of function evaluations=2 >>> KSP Object: 1 MPI processes >>> type: gmres >>> GMRES: restart=30, using Classical (unmodified) GramSchmidt >>> Orthogonalization with no iterative refinement >>> GMRES: happy breakdown tolerance 1e30 >>> maximum iterations=10000, initial guess is zero >>> tolerances: relative=1e05, absolute=1e50, divergence=10000 >>> left preconditioning >>> using PRECONDITIONED norm type for convergence test >>> PC Object: 1 MPI processes >>> type: ilu >>> ILU: outofplace factorization >>> 0 levels of fill >>> tolerance for zero pivot 2.22045e14 >>> using diagonal shift to prevent zero pivot >>> matrix ordering: natural >>> factor fill ratio given 1, needed 1 >>> Factored matrix follows: >>> Matrix Object: 1 MPI processes >>> type: seqaij >>> rows=153, cols=153 >>> package used to perform factorization: petsc >>> total: nonzeros=2401, allocated nonzeros=2401 >>> total number of mallocs used during MatSetValues calls =0 >>> not using Inode routines >>> linear system matrix = precond matrix: >>> Matrix Object: 1 MPI processes >>> type: seqaij >>> rows=153, cols=153 >>> total: nonzeros=2401, allocated nonzeros=2401 >>> total number of mallocs used during MatSetValues calls =0 >>> not using Inode routines >>> SNESLineSearch Object: 1 MPI processes >>> type: basic >>> maxstep=1e+08, minlambda=1e12 >>> tolerances: relative=1e08, absolute=1e15, lambda=1e08 >>> maximum iterations=1 >>> Fracture system solved at nonlinear iteration 1 , final nonlinear >>> residual norm: 5.10144e07 >>> >>>  without the flag  >>> >>> applying variable bounds ... done >>> NL step 0, residual_2 = 2.41046e01 >>> SNES Jacobian assembly... done >>> Linear solve converged due to CONVERGED_RTOL iterations 3 >>> NL step 1, residual_2 = 5.10144e07 >>> SNES Jacobian assembly... done >>> Linear solve converged due to CONVERGED_RTOL iterations 3 >>> NL step 2, residual_2 = 1.61153e12 >>> SNES Object: 1 MPI processes >>> type: virs >>> maximum iterations=50, maximum function evaluations=10000 >>> tolerances: relative=1e08, absolute=1e35, solution=1e08 >>> total number of linear solver iterations=6 >>> total number of function evaluations=3 >>> KSP Object: 1 MPI processes >>> type: gmres >>> GMRES: restart=30, using Classical (unmodified) GramSchmidt >>> Orthogonalization with no iterative refinement >>> GMRES: happy breakdown tolerance 1e30 >>> maximum iterations=10000, initial guess is zero >>> tolerances: relative=1e05, absolute=1e50, divergence=10000 >>> left preconditioning >>> using PRECONDITIONED norm type for convergence test >>> PC Object: 1 MPI processes >>> type: ilu >>> ILU: outofplace factorization >>> 0 levels of fill >>> tolerance for zero pivot 2.22045e14 >>> using diagonal shift to prevent zero pivot >>> matrix ordering: natural >>> factor fill ratio given 1, needed 1 >>> Factored matrix follows: >>> Matrix Object: 1 MPI processes >>> type: seqaij >>> rows=153, cols=153 >>> package used to perform factorization: petsc >>> total: nonzeros=2401, allocated nonzeros=2401 >>> total number of mallocs used during MatSetValues calls =0 >>> not using Inode routines >>> linear system matrix = precond matrix: >>> Matrix Object: 1 MPI processes >>> type: seqaij >>> rows=153, cols=153 >>> total: nonzeros=2401, allocated nonzeros=2401 >>> total number of mallocs used during MatSetValues calls =0 >>> not using Inode routines >>> SNESLineSearch Object: 1 MPI processes >>> type: basic >>> maxstep=1e+08, minlambda=1e12 >>> tolerances: relative=1e08, absolute=1e15, lambda=1e08 >>> maximum iterations=1 >>> Fracture system solved at nonlinear iteration 2 , final nonlinear >>> residual norm: 1.61153e12 >>> >>> >>> >>> Incidentally, there _was_ a reference counting bug that led to a >>> breakdown very much like what you reported. >>> I'm attaching that patch and a small improvement to DMVariableBounds. I >>> need to figure out how to go about creating matrices from DM_libMesh. >>> >>> Let me know if you found anything or if I can be of any help. >>> >>> Thanks, >>> Ata >>> >>> Dmitry. >>> >>> On Tue, May 1, 2012 at 6:21 PM, Ataollah Mesgarnejad < >>> amesga1@...> wrote: >>> >>>> Dmitry, >>>> >>>> Is there a way so SNES doesn't go calculate the Jacobian for every >>>> iteration (my Jac does't change)? I used to use a flag to check if it's >>>> assembled but when I do that with petscdm, KSP associated with SNES gives >>>> DIVERGED_BREAKDOWN for the second iteration! looks like petscdm doesn't >>>> make a copy of the matrix so when it's used in KSP solve it gets jumbled up. >>>> >>>> Here is how it looks >>>> >>>> with the flag: >>>> >>>> applying variable bounds ... done >>>> NL step 0, residual_2 = 8.90674e04 >>>> SNES Jacobian assembly... done >>>> Linear solve converged due to CONVERGED_RTOL iterations 6 >>>> NL step 1, residual_2 = 2.90987e05 >>>> Linear solve did not converge due to DIVERGED_BREAKDOWN iterations 1 >>>> >>>> this happens for every SNES solve. >>>>  >>>> >>>> without it: >>>> >>>> applying variable bounds ... done >>>> NL step 0, residual_2 = 6.08030e02 >>>> SNES hessian assembly... done >>>> Linear solve converged due to CONVERGED_RTOL iterations 6 >>>> NL step 1, residual_2 = 2.61557e07 >>>> SNES hessian assembly... done >>>> Linear solve converged due to CONVERGED_RTOL iterations 6 >>>> NL step 2, residual_2 = 1.35354e12 >>>> Fracture system solved at nonlinear iteration 2 , final nonlinear >>>> residual norm: 1.35354e12 >>>> >>>> ps: these are not from the same time step but you get the idea. >>>> >>>> Best, >>>> Ata >>>> >>>> On May 1, 2012, at 8:49 AM, Dmitry Karpeev wrote: >>>> >>>> As an fyi: I just built and ran miscellaneous_ex7 in parallel with the >>>> latest petscdev: >>>> http://petsc.cs.iit.edu/petsc/petscdev/rev/12d4a1ea0aca >>>> If you decide you want to upgrade to this (or later) petscdev, you'd >>>> need a small patch for >>>> petsc_dm_nonlinear_solver.C to account for the latest API >>>> changes (attached). >>>> I am working on a bigger patch that adds more functionality to >>>> DM_libMesh and streamline >>>> miscellaneous_ex7. That will take some time to percolate up into the >>>> libMesh tree. In the >>>> meantime you should be able to run with the small fix I've attached. >>>> >>>> Let me know how it works. >>>> Dmitry. >>>> >>>> >>>> On Tue, May 1, 2012 at 8:35 AM, Ataollah Mesgarnejad < >>>> amesga1@...> wrote: >>>> >>>>> Dmitry, >>>>> >>>>> I believe that the error is indeed due to petscdev version. I'm going >>>>> to check this during the next few hours and get back to you. >>>>> >>>>> Thanks, >>>>> Ata >>>>> >>>>> >>>>> On Tue, May 1, 2012 at 8:29 AM, Dmitry Karpeev <karpeev@...>wrote: >>>>> >>>>>> Ata, >>>>>> I can't reproduce the problem: miscellaneous_ex7 runs fine for me on >>>>>> 2 procs with a relatively recent version >>>>>> of petscdev. Can you send me the full PETSc error output? >>>>>> >>>>>> Thanks. >>>>>> Dmitry. >>>>>> >>>>>> >>>>>> On Mon, Apr 30, 2012 at 8:56 PM, Ataollah Mesgarnejad < >>>>>> amesga1@...> wrote: >>>>>> >>>>>>> >>>>>>> On Apr 30, 2012, at 8:39 PM, Dmitry Karpeev wrote: >>>>>>> >>>>>>> Ata, >>>>>>> example7 segfaults when you run in parallel? >>>>>>> Does it do it on 2 procs? >>>>>>> >>>>>>> Yes and yes. >>>>>>> >>>>>>> here is the stack: >>>>>>> (gdb) where >>>>>>> #0 0x00000038fbaaada0 in __nanosleep_nocancel () from >>>>>>> /lib64/libc.so.6 >>>>>>> #1 0x00000038fbaaac30 in sleep () from /lib64/libc.so.6 >>>>>>> #2 0x00002ae457b35c40 in PetscSleep(double) () >>>>>>> from /share/apps/petscdev/Linuxintel12.1mef90O/lib/libpetsc.so >>>>>>> #3 0x00002ae457b7dffa in PetscAttachDebugger() () >>>>>>> from /share/apps/petscdev/Linuxintel12.1mef90O/lib/libpetsc.so >>>>>>> #4 0x00002ae457b7e95f in PetscAttachDebuggerErrorHandler(int, int, >>>>>>> char const*, char const*, char const*, int, PetscErrorType, char const*, >>>>>>> void*) () >>>>>>> from /share/apps/petscdev/Linuxintel12.1mef90O/lib/libpetsc.so >>>>>>> #5 0x00002ae457b817be in PetscError(int, int, char const*, char >>>>>>> const*, char const*, int, PetscErrorType, char const*, ...) () >>>>>>> from /share/apps/petscdev/Linuxintel12.1mef90O/lib/libpetsc.so >>>>>>> #6 0x00002ae457b82ea7 in PetscDefaultSignalHandler(int, void*) () >>>>>>> from /share/apps/petscdev/Linuxintel12.1mef90O/lib/libpetsc.so >>>>>>> #7 0x00002ae457b82883 in PetscSignalHandler_Private () >>>>>>> from /share/apps/petscdev/Linuxintel12.1mef90O/lib/libpetsc.so >>>>>>> #8 <signal handler called> >>>>>>> #9 0x00002ae458189340 in DMCoarsenHookAdd(_p_DM*, int (*)(_p_DM*, >>>>>>> _p_DM*, void*), int (*)(_p_DM*, _p_Mat*, _p_Vec*, _p_Mat*, _p_DM*, void*), >>>>>>> void*) () >>>>>>> from /share/apps/petscdev/Linuxintel12.1mef90O/lib/libpetsc.so >>>>>>> #10 0x00002ae45835f431 in DMSNESGetContext(_p_DM*, _n_SNESDM**) () >>>>>>> from /share/apps/petscdev/Linuxintel12.1mef90O/lib/libpetsc.so >>>>>>> #11 0x00002ae45835ea34 in DMSNESGetFunction(_p_DM*, int >>>>>>> (**)(_p_SNES*, _p_Vec*, Type <return> to continue, or q <return> to >>>>>>> quit >>>>>>> _p_Vec*, void*), void**) () >>>>>>> from /share/apps/petscdev/Linuxintel12.1mef90O/lib/libpetsc.so >>>>>>> #12 0x00002ae458393d1c in SNESGetFunction(_p_SNES*, _p_Vec**, int >>>>>>> (**)(_p_SNES*, _p_Vec*, _p_Vec*, void*), void**) () >>>>>>> from /share/apps/petscdev/Linuxintel12.1mef90O/lib/libpetsc.so >>>>>>> #13 0x00002ae45839b9eb in SNESSetUp(_p_SNES*) () >>>>>>> from /share/apps/petscdev/Linuxintel12.1mef90O/lib/libpetsc.so >>>>>>> #14 0x00002ae45839a6b0 in SNESSolve(_p_SNES*, _p_Vec*, _p_Vec*) () >>>>>>> from /share/apps/petscdev/Linuxintel12.1mef90O/lib/libpetsc.so >>>>>>> >>>>>>> #15 0x00000000010976ab in >>>>>>> libMesh::PetscDMNonlinearSolver<libMesh::Number>::solve (this=0x35506d0, >>>>>>> jac_in=..., x_in=..., r_in=...) >>>>>>> at src/solvers/petsc_dm_nonlinear_solver.C:450 >>>>>>> #16 0x0000000000cd604f in libMesh::NonlinearImplicitSystem::solve ( >>>>>>> this=0x353b7e0) at src/systems/nonlinear_implicit_system.C:178 >>>>>>> #17 0x0000000000426082 in Biharmonic::step (this=0x34f2000, >>>>>>> dt=@0x7fff95edca28) >>>>>>> at biharmonic.C:110 >>>>>>> #18 0x0000000000426373 in Biharmonic::run (this=0x34f2000) at >>>>>>> biharmonic.C:152 >>>>>>> warning: Range for type (null) has invalid bounds 0..111 >>>>>>> #19 0x0000000000436cf5 in main (argc=12, argv=0x7fff95edcc18) >>>>>>> at miscellaneous_ex7.C:64 >>>>>>> >>>>>>> Ata >>>>>>> >>>>>>> >>>>>>> Dmitry. >>>>>>> >>>>>>> On Mon, Apr 30, 2012 at 8:31 PM, Ataollah Mesgarnejad < >>>>>>> amesga1@...> wrote: >>>>>>> >>>>>>>> Dmitry, >>>>>>>> >>>>>>>> I had a very good experience on single processor but on multiple >>>>>>>> processors I get segmentation error at line where we call solve (I tried >>>>>>>> both my code and example7) . I tried your example without vi >>>>>>>> usepetscdm and it works just fine. Right now I'm trying to compile a >>>>>>>> debug version and look into it; meanwhile I thought it would be a good idea >>>>>>>> if I let you know. >>>>>>>> >>>>>>>> Thanks, >>>>>>>> Ata >>>>>>>> >>>>>>>> On Apr 30, 2012, at 1:30 PM, Dmitry Karpeev wrote: >>>>>>>> >>>>>>>> >>>>>>>> >>>>>>>> On Mon, Apr 30, 2012 at 1:28 PM, Ataollah Mesgarnejad < >>>>>>>> amesga1@...> wrote: >>>>>>>> >>>>>>>>> >>>>>>>>> On Apr 30, 2012, at 1:13 PM, Dmitry Karpeev wrote: >>>>>>>>> >>>>>>>>> >>>>>>>>> >>>>>>>>> On Mon, Apr 30, 2012 at 1:03 PM, Ataollah Mesgarnejad < >>>>>>>>> amesga1@...> wrote: >>>>>>>>> >>>>>>>>>> Dmitry, >>>>>>>>>> >>>>>>>>>> I had time to come back and work on the SNESVI/Libmesh program. >>>>>>>>>> The way you setup your code is particularly hard for me to follow in my own >>>>>>>>>> code since it requires a significant amount of extra work (to define >>>>>>>>>> classes and setup the inheritances..). What I did instead was to try add a >>>>>>>>>> TransientNonlinearSystem to my EquationSystems object: >>>>>>>>>> >>>>>>>>>> TransientNonlinearImplicitSystem & fracture_system = >>>>>>>>>> equation_system.add_system<TransientNonlinearImplicitSystem>("fracture_system"); >>>>>>>>>> >>>>>>>>>> and then to add the bounds and Objective and Residual functions: >>>>>>>>>> >>>>>>>>>> fracture_system.attach_init_function(init_ic_fracture); >>>>>>>>>> fracture_system.nonlinear_solver>jacobian= PsiFormJacobian; >>>>>>>>>> fracture_system.nonlinear_solver>residual= PsiFormResidual; >>>>>>>>>> fracture_system.nonlinear_solver>bounds = PsiFormBounds; >>>>>>>>>> >>>>>>>>>> what I see now is that SNES solve does not apply the bounds. So >>>>>>>>>> here is my questions: >>>>>>>>>> >>>>>>>>>> 1 Is there inherently wrong with setting up the system >>>>>>>>>> specially the bounds the easy way like I did? >>>>>>>>>> >>>>>>>>> Ata, >>>>>>>>> Doing things this "simple" way is just fine: there is no reason to >>>>>>>>> follow the relatively convoluted structure I used in that >>>>>>>>> example. When I have time, I'll rewrite it to simplify things. >>>>>>>>> >>>>>>>>>> 2 Where is the right SNES type is chosen ? I guess this should >>>>>>>>>> be hard wired somewhere in NonlinearImplicitSystem! >>>>>>>>>> 3 When I use snes_view (with your example or my own code ) I >>>>>>>>>> get snes_type ls? shouldn't this be VI ? >>>>>>>>>> >>>>>>>>> >>>>>>>>> What's likely missing is the commandline option usepetscdm. >>>>>>>>> It's necessary since libMesh wants to preserve the usual behavior as >>>>>>>>> the default. In particular, because SNESVI and the DM hookup aren't really >>>>>>>>> possible prior to petsc3.2. >>>>>>>>> In fact, it might not be possible without a recent petscdev (let >>>>>>>>> me verify that, once I get to my computer). >>>>>>>>> >>>>>>>>> In the meantime, could you give usepetscdm a try? >>>>>>>>> >>>>>>>>> I just saw that in the example help line and used it. It works >>>>>>>>> like a charm. >>>>>>>>> >>>>>>>> Great! >>>>>>>> Dmitry. >>>>>>>> >>>>>>>>> I'll go on check for the convergence and see what I can get. >>>>>>>>> >>>>>>>>> Thanks, >>>>>>>>> Ata >>>>>>>>> >>>>>>>>> Thanks. >>>>>>>>> >>>>>>>>> Dmitry. >>>>>>>>> >>>>>>>>>> >>>>>>>>>> Thanks, >>>>>>>>>> Ata >>>>>>>>>> >>>>>>>>>> On Apr 12, 2012, at 10:00 AM, Dmitry Karpeev wrote: >>>>>>>>>> >>>>>>>>>> Convergence is something we might need to work on  I noticed >>>>>>>>>> some chatter in a simple bounded biharmonic heat equation (a stripped down >>>>>>>>>> version of CahnHilliard). >>>>>>>>>> Thanks. >>>>>>>>>> Dmitry. >>>>>>>>>> >>>>>>>>>> On Thu, Apr 12, 2012 at 9:55 AM, Ataollah Mesgarnejad < >>>>>>>>>> amesga1@...> wrote: >>>>>>>>>> >>>>>>>>>>> Thanks Dmitry, >>>>>>>>>>> >>>>>>>>>>> It would be wonderful. I made a wrapper for SNESVI in myself for >>>>>>>>>>> my program but had difficulty achieving good convergence. Right now I am >>>>>>>>>>> occupied with other issues but I'll get back to this as soon as I get a >>>>>>>>>>> chance. >>>>>>>>>>> >>>>>>>>>>> Ata >>>>>>>>>>> >>>>>>>>>>> On Apr 10, 2012, at 9:55 PM, Dmitry Karpeev wrote: >>>>>>>>>>> >>>>>>>>>>> Ata, >>>>>>>>>>> There is now a way to use SNESVI from libMesh. >>>>>>>>>>> If interested, pull a recent revision from sourceforge and take >>>>>>>>>>> a look at examples/miscellaneous/miscellaneous_ex7. >>>>>>>>>>> Let me know, if you have any questions. >>>>>>>>>>> >>>>>>>>>>> Thanks. >>>>>>>>>>> Dmitry. >>>>>>>>>>> >>>>>>>>>>> On Wed, Dec 14, 2011 at 12:45 PM, Dmitry Karpeev < >>>>>>>>>>> karpeev@...> wrote: >>>>>>>>>>> >>>>>>>>>>>> Okay, thanks for your feedback. >>>>>>>>>>>> I'll keep you updated on this work. >>>>>>>>>>>> >>>>>>>>>>>> Dmitry. >>>>>>>>>>>> >>>>>>>>>>>> >>>>>>>>>>>> On Wed, Dec 14, 2011 at 12:42 PM, Ataollah Mesgarnejad < >>>>>>>>>>>> amesga1@...> wrote: >>>>>>>>>>>> >>>>>>>>>>>>> >>>>>>>>>>>>> >>>>>>>>>>>>> On Wed, Dec 14, 2011 at 12:35 PM, Dmitry Karpeev < >>>>>>>>>>>>> karpeev@...> wrote: >>>>>>>>>>>>> >>>>>>>>>>>>>> >>>>>>>>>>>>>> >>>>>>>>>>>>>> On Wed, Dec 14, 2011 at 12:19 PM, Ataollah Mesgarnejad < >>>>>>>>>>>>>> amesga1@...> wrote: >>>>>>>>>>>>>> >>>>>>>>>>>>>>> >>>>>>>>>>>>>>> >>>>>>>>>>>>>>> On Wed, Dec 14, 2011 at 11:59 AM, Dmitry Karpeev < >>>>>>>>>>>>>>> karpeev@...> wrote: >>>>>>>>>>>>>>> >>>>>>>>>>>>>>>> Ata, >>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>> A couple of things: >>>>>>>>>>>>>>>> 1. SNESVI solves a VI under what's known as "box >>>>>>>>>>>>>>>> constraints": the Vec x that the residual operates on is entrywise between >>>>>>>>>>>>>>>> Vecs xlow and xhigh. In order to enable >>>>>>>>>>>>>>>> constraints on a linear combination of the Vec entries (as >>>>>>>>>>>>>>>> you suggest), the problem >>>>>>>>>>>>>>>> has to be augmented with extra variables  these linear >>>>>>>>>>>>>>>> combinations. This, of course, >>>>>>>>>>>>>>>> changes the Jacobian, preconditioner, etc. The user >>>>>>>>>>>>>>>> provides only the Jacobian for the original system, so libMesh would have >>>>>>>>>>>>>>>> to automatically construct the Jacobian for the extended system, and >>>>>>>>>>>>>>>> precondition it using the user provided PC for the original subsystem. This >>>>>>>>>>>>>>>> is not necessarily the unreasonable thing to do in many circumstances, but >>>>>>>>>>>>>>>> will take more thinking than a straighforward nodal constraint. The user >>>>>>>>>>>>>>>> API for specifying the constraints, however, becomes rather straighforward. >>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>> 2. In terms of FEM convergence, constraining at the nodes >>>>>>>>>>>>>>>> is sufficient (see, for example, >>>>>>>>>>>>>>>> Glowinski's book: >>>>>>>>>>>>>>>> http://books.google.com/books/about/Numerical_analysis_of_variational_inequa.html?id=Pf4ed2mtbx4C >>>>>>>>>>>>>>>> ) >>>>>>>>>>>>>>>> You can show that if you impose the nodal constraints, in >>>>>>>>>>>>>>>> the limit of mesh size going to 0 the inequality is satisfied a.e. >>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>> In this case no funny manipulation of the Jacobian or the >>>>>>>>>>>>>>>> PC is necessary, but the user API to specify the constraint is a bit >>>>>>>>>>>>>>>> amiguous: how do you constrain the node? From a containing element? Then >>>>>>>>>>>>>>>> there are multiple, possibly incompatible, constraints. >>>>>>>>>>>>>>>> >>>>>>>>>>>>>>> Or you can simply go through all the nodes in the system >>>>>>>>>>>>>>> once using node iterator. Here is how I think you can do it: >>>>>>>>>>>>>>> >>>>>>>>>>>>>>>  >>>>>>>>>>>>>>> //make two PetscVector ub,lb from two variables you already >>>>>>>>>>>>>>> added to the system >>>>>>>>>>>>>>> NumericVector<Number> &ub_sys = >>>>>>>>>>>>>>> fracture_system.get_vector("ub"); >>>>>>>>>>>>>>> NumericVector<Number> &lb_sys = >>>>>>>>>>>>>>> fracture_system.get_vector("lb"); >>>>>>>>>>>>>>> PetscVector<Number>&lb = >>>>>>>>>>>>>>> libmesh_cast_ptr<PetscVector<Number>*>(&ub_sys); >>>>>>>>>>>>>>> PetscVector<Number>&ub = >>>>>>>>>>>>>>> libmesh_cast_ptr<PetscVector<Number>*>(&lb_sys); >>>>>>>>>>>>>>> >>>>>>>>>>>>>>> *then in the bound routine do:* >>>>>>>>>>>>>>> >>>>>>>>>>>>>>> // Get a constant reference to the mesh object. >>>>>>>>>>>>>>> const MeshBase& mesh = es.get_mesh(); >>>>>>>>>>>>>>> //Get the iterator >>>>>>>>>>>>>>> MeshBase::const_node_iterator nd = >>>>>>>>>>>>>>> mesh.local_nodes_begin(); >>>>>>>>>>>>>>> const MeshBase::const_node_iterator end_nd = >>>>>>>>>>>>>>> mesh.local_nodes_end(); >>>>>>>>>>>>>>> for (; nd != end_nd; nd++){ >>>>>>>>>>>>>>> const Node* node = *nd; >>>>>>>>>>>>>>> //get the degree of freedom of node >>>>>>>>>>>>>>> unsigned int dn = node>dof_number(0,0,0); // >>>>>>>>>>>>>>> system 0 , var 0, component 0 >>>>>>>>>>>>>>> //set the ub,lb >>>>>>>>>>>>>>> ub(dn)= A; >>>>>>>>>>>>>>> lb(dn) = B; >>>>>>>>>>>>>>> } >>>>>>>>>>>>>>> //close the vectors >>>>>>>>>>>>>>> ub.close(); >>>>>>>>>>>>>>> lb.close(); >>>>>>>>>>>>>>> >>>>>>>>>>>>>> >>>>>>>>>>>>>> You could do that. That's probably the mode we can support >>>>>>>>>>>>>> almost immediately. >>>>>>>>>>>>>> The issue here is that you lose some of the "context" by >>>>>>>>>>>>>> processing each node by >>>>>>>>>>>>>> itself rather as part of an element. For example what if >>>>>>>>>>>>>> each element is a different material? >>>>>>>>>>>>>> I guess in that case that property should be set on the nodes >>>>>>>>>>>>>> as well. >>>>>>>>>>>>>> >>>>>>>>>>>>> >>>>>>>>>>>>> Exactly one can always pass mat prop as a Constant MONOMIAL >>>>>>>>>>>>> (e.g. system.add_variable("mat_prop",CONSTANT,MONOMIAL);) variable with >>>>>>>>>>>>> values on each node. >>>>>>>>>>>>> >>>>>>>>>>>>> Ata >>>>>>>>>>>>> >>>>>>>>>>>>> >>>>>>>>>>>>>> Dmitry. >>>>>>>>>>>>>> >>>>>>>>>>>>>>> >>>>>>>>>>>>>>>  >>>>>>>>>>>>>>> Ata >>>>>>>>>>>>>>> >>>>>>>>>>>>>>> >>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>> Dmitry. >>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>> On Wed, Dec 14, 2011 at 11:17 AM, Ataollah Mesgarnejad < >>>>>>>>>>>>>>>> amesga1@...> wrote: >>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>> I hope I understood the whole thing correctly. >>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>> On Wed, Dec 14, 2011 at 10:56 AM, Dmitry Karpeev < >>>>>>>>>>>>>>>>> karpeev@...> wrote: >>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>> On Wed, Dec 14, 2011 at 10:50 AM, Ataollah Mesgarnejad < >>>>>>>>>>>>>>>>>> amesga1@...> wrote: >>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>> I don't think setting the variable bound on each >>>>>>>>>>>>>>>>>>> node separately is a good idea in general. >>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>> I'm not sure what you mean by this. Could you clarify? >>>>>>>>>>>>>>>>>> I don't think I'm advocating this. >>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>> Bounds should be set as a sum on each elements degrees of >>>>>>>>>>>>>>>>>>> freedom; this way there should not be any nodal inconsistencies. >>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>> I'm not sure what you mean here either. A bound is >>>>>>>>>>>>>>>>>> usually not a sum (e.g., a bound on a concentration is 0 <= c <=1). >>>>>>>>>>>>>>>>>> Perhaps your application expresses the total bound as a >>>>>>>>>>>>>>>>>> sum? I'd be very interested in learning more about that case. >>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>> OK here is my thought process: one optimally wants to >>>>>>>>>>>>>>>>> bound the variable over the whole Element domain between lb,ub. So we would >>>>>>>>>>>>>>>>> have something like this: >>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>> lb \leq \sum_{i=1}^{n=NDOF} \phi_{i}(x) \hat{u}_{i} \leq >>>>>>>>>>>>>>>>> ub, \forall x \in El >>>>>>>>>>>>>>>>> (Eq1) >>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>> so it's not exactly that the nodal values should be >>>>>>>>>>>>>>>>> bounded (\hat{u}_{i}) unless it's a FD application where \phi_i(x) is delta >>>>>>>>>>>>>>>>> function !! That said I think given a normalized basis, bounding nodal >>>>>>>>>>>>>>>>> values is stricter requirement that (Eq1). if that's what you meant by >>>>>>>>>>>>>>>>> inconsistencies between node bounds then the stronger bound would be my >>>>>>>>>>>>>>>>> first priority. >>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>> Ata >>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>> Of course it might run into a situations with >>>>>>>>>>>>>>>>>>> semismooth VI. >>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>> The type of VI solver is independent of how the bounds >>>>>>>>>>>>>>>>>> are specified: once we know the upper and lower bounds for each element of >>>>>>>>>>>>>>>>>> a Vec, we can solve the VI with either a semismooth method or an active set >>>>>>>>>>>>>>>>>> method. >>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>> I really appreciate your feedback. I hope we can get to >>>>>>>>>>>>>>>>>> the point where things are as clear as possible so that you can >>>>>>>>>>>>>>>>>> profitably use the new capabilities. >>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>> Dmitry. >>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>> Thanks, >>>>>>>>>>>>>>>>>>> Ata >>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>> On Wed, Dec 14, 2011 at 9:47 AM, Dmitry Karpeev < >>>>>>>>>>>>>>>>>>> karpeev@...> wrote: >>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>>> Are you familiar with SNESVI? The main thing I'm >>>>>>>>>>>>>>>>>>>> looking at right now is the libMesh API for specifying >>>>>>>>>>>>>>>>>>>> the bounds on the variables. It seems to me that the >>>>>>>>>>>>>>>>>>>> user would ideally want to specify the bounds on the variable >>>>>>>>>>>>>>>>>>>> one element at a time, since that's the default >>>>>>>>>>>>>>>>>>>> traversal mode for libMesh. Then, of course, many nodal values >>>>>>>>>>>>>>>>>>>> are traversed many time, allowing for an inconsistency. >>>>>>>>>>>>>>>>>>>> The question is, what you, as a user, expect in that situation? >>>>>>>>>>>>>>>>>>>> Should the latest bound that the user has set be >>>>>>>>>>>>>>>>>>>> applied or the tightest of the bounds the user has supplied? >>>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>>> Dmitry. >>>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>>> On Wed, Dec 14, 2011 at 9:28 AM, Ataollah Mesgarnejad < >>>>>>>>>>>>>>>>>>>> amesga1@...> wrote: >>>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>>>> Thanks Dmitry, I look forward to use it. Also if it's >>>>>>>>>>>>>>>>>>>>> available as a dev version I would love to try to work with it. >>>>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>>>> Ata >>>>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>>>> On Wed, Dec 14, 2011 at 8:24 AM, Dmitry Karpeev < >>>>>>>>>>>>>>>>>>>>> karpeev@...> wrote: >>>>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>>>>> We (ANL) are working on an patch to libMesh to enable >>>>>>>>>>>>>>>>>>>>>> this functionality. >>>>>>>>>>>>>>>>>>>>>> The motivation is to use it through MOOSE, but it >>>>>>>>>>>>>>>>>>>>>> will be exposed as a libMesh API. >>>>>>>>>>>>>>>>>>>>>> Dmitry. >>>>>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>>>>> On Wed, Dec 14, 2011 at 7:53 AM, Ataollah Mesgarnejad >>>>>>>>>>>>>>>>>>>>>> <amesga1@...> wrote: >>>>>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>>>>>> Dear all, >>>>>>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>>>>>> I wonder if there is compatible versionIs it >>>>>>>>>>>>>>>>>>>>>>> possible to >>>>>>>>>>>>>>>>>>>>>>> use PETScNonlinearSolver class with SNESVI >>>>>>>>>>>>>>>>>>>>>>> (PETSc/3.2) to solve a >>>>>>>>>>>>>>>>>>>>>>> variational inequality problem? >>>>>>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>>>>>> Thanks, >>>>>>>>>>>>>>>>>>>>>>> Ata Mesgarnejad >>>>>>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>>>>>>  >>>>>>>>>>>>>>>>>>>>>>> Cloud Computing  Latest Buzzword or a Glimpse of >>>>>>>>>>>>>>>>>>>>>>> the Future? >>>>>>>>>>>>>>>>>>>>>>> This paper surveys cloud computing today: What are >>>>>>>>>>>>>>>>>>>>>>> the benefits? >>>>>>>>>>>>>>>>>>>>>>> Why are businesses embracing it? What are its >>>>>>>>>>>>>>>>>>>>>>> payoffs and pitfalls? >>>>>>>>>>>>>>>>>>>>>>> http://www.accelacomm.com/jaw/sdnl/114/51425149/ >>>>>>>>>>>>>>>>>>>>>>> _______________________________________________ >>>>>>>>>>>>>>>>>>>>>>> Libmeshusers mailing list >>>>>>>>>>>>>>>>>>>>>>> Libmeshusers@... >>>>>>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>>>>>> https://lists.sourceforge.net/lists/listinfo/libmeshusers >>>>>>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>>>>  >>>>>>>>>>>>>>>>>>>>> A. Mesgarnejad >>>>>>>>>>>>>>>>>>>>> PhD Student, Research Assistant >>>>>>>>>>>>>>>>>>>>> Mechanical Engineering Department >>>>>>>>>>>>>>>>>>>>> Louisiana State University >>>>>>>>>>>>>>>>>>>>> 2203 Patrick F. Taylor Hall >>>>>>>>>>>>>>>>>>>>> Baton Rouge, La 70803 >>>>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>>  >>>>>>>>>>>>>>>>>>> A. Mesgarnejad >>>>>>>>>>>>>>>>>>> PhD Student, Research Assistant >>>>>>>>>>>>>>>>>>> Mechanical Engineering Department >>>>>>>>>>>>>>>>>>> Louisiana State University >>>>>>>>>>>>>>>>>>> 2203 Patrick F. Taylor Hall >>>>>>>>>>>>>>>>>>> Baton Rouge, La 70803 >>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>> >>>>>>>>>>>>>>> >>>>>>>>>>>>>>> >>>>>>>>>>>>>>>  >>>>>>>>>>>>>>> A. Mesgarnejad >>>>>>>>>>>>>>> PhD Student, Research Assistant >>>>>>>>>>>>>>> Mechanical Engineering Department >>>>>>>>>>>>>>> Louisiana State University >>>>>>>>>>>>>>> 2203 Patrick F. Taylor Hall >>>>>>>>>>>>>>> Baton Rouge, La 70803 >>>>>>>>>>>>>>> >>>>>>>>>>>>>> >>>>>>>>>>>>>> >>>>>>>>>>>>> >>>>>>>>>>>>> >>>>>>>>>>>>>  >>>>>>>>>>>>> A. Mesgarnejad >>>>>>>>>>>>> PhD Student, Research Assistant >>>>>>>>>>>>> Mechanical Engineering Department >>>>>>>>>>>>> Louisiana State University >>>>>>>>>>>>> 2203 Patrick F. Taylor Hall >>>>>>>>>>>>> Baton Rouge, La 70803 >>>>>>>>>>>>> >>>>>>>>>>>> >>>>>>>>>>>> >>>>>>>>>>> >>>>>>>>>>> >>>>>>>>>> >>>>>>>>>> >>>>>>>>> >>>>>>>>> >>>>>>>> >>>>>>>> >>>>>>> >>>>>>> >>>>>> >>>>> >>>>> >>>>>  >>>>> A. Mesgarnejad >>>>> PhD Student, Research Assistant >>>>> Mechanical Engineering Department >>>>> Louisiana State University >>>>> 2203 Patrick F. Taylor Hall >>>>> Baton Rouge, La 70803 >>>>> >>>> >>>> <petsc_dm.patch> >>>> >>>> >>>> >>> <0001PETScDMincrefbugfix.patch> >>> <0002PETScDMinitializesboundswithInfnow.patch> >>> >>> >>> >> >> > > 