From: Dmitry Karpeev <karpeev@mc...>  20120505 17:22:06

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