[petsc-users] Conditional Constraints

[Shri]

Jon,
    I am trying to compile petsc with complex data type right now and will try to test the VI solver with complex variables within the next few hours. 


> My understanding is that comparison operators (such as < and >) are
> not defined for complex numbers, so would it be correct to say that
> anyone using complex scalars within PETSc will encounter this problem
> with SNESVI? In other words, SNESVI inequality constraints can only be
> applied to real variables or real-valued functions of complex
> variables?

For complex scalars, we only compare the real part of the complex variable and bound,i.e.,
real(xl) <= real(x) <= real(x_u). For your application, this real part comparison is sufficient 
since the magnitude is saved in the real part.


> >> The residuals are set as follows, and the Jacobian entries are set
> >> correspondingly:
> >> f[V_block] = {function of adjacent V_blocks, V_electrode}
> >> f[V_electrode] = 0.0 + 0.0i
> >> f[magnitude(V_electrode)] = 0.0 + 0.0i
> >> f[magnitude(I_electrode)] = 0.0 + 0.0i

Can you give an example of f for all dofs maybe just for one grid point with the following dof names
Vb - V_block
Ve - V_electrode
Vmag - magnitude(V_electrode)
Imag - magnitude(I_electrode) 

> Since I set all the residuals for both magnitude variables to zero,
> and each entry in the Jacobian matrix is a partial derivative of a
> residual, I believe the Jacobian entries for both magnitudes would all
> be zero as well. However, your point made me realize I have no
> Jacobian entry for the partial derivative of f[V_block] w.r.t.
> V_electrode. Is that needed? 

Since f[V_block] is a function of V_electrode, you would need to set the partial derivatives of f[V_block] w.r.t V_electrode
for correct jacobian evaluation and for having a good newton convergence.

How would one enter it with
> MatSetValuesStencil? 
ex28.c,
http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-dev/src/ksp/ksp/examples/tutorials/ex28.c.html,
is a petsc example which shows how to use MatSetValuesStencil with multiple degrees of freedom. 
See the routine ComputeMatrix in the example.

My working understanding was that each degree of
> freedom conceptually has its own residual and Jacobian, even though
> all residuals are contained in one Vec and all Jacobians contained in
> one Mat.
> 

Please bare with us as we are trying to figure out whether the issue is with our VI implementation or something in 
your code is amiss. I appreciate your patience on this matter.

Btw, one quick question, are all the problem formulations in your field done using complex variables or by decomposing the complex term into real/imaginary or magnitude/angle? My background is in electrical power grid and we mostly use either real/imaginary part or magnitude/angle as the variables.


Shri


----- Original Message -----
> Thanks for your reply, Shri. Responses below:
> 
> On 2011-07-28, at 1:03 PM, Shri wrote:
> 
> >
> > ----- Original Message -----
> >> Hi Barry, Shri,
> >>
> >> I reworked my implementation to include only the time-independent
> >> electrical problem, and to make use of SNESVISetVariableBounds().
> >> However, it does not work quite as I expected: Constraints on
> >> variables that are functions of other variables have no effect on
> >> the
> >> variables by which the functions are defined.
> >>
> >> My constraints are set on the four degrees of freedom as follows:
> >> -SNES_VI_INF <= V_block <= SNES_VI_INF
> >> -SNES_VI_INF <= V_electrode <= SNES_VI_INF
> >> 0.0 + 0.0i <= magnitude(V_electrode) <= magnitude(V_electrode_max)
> >> 0.0 + 0.0i <= magnitude(I_electrode) <= magnitude(I_electrode_max)
> >
> > Note that while V_block and V_electrode and complex quantities,the
> > magnitude terms are actually real numbers. What function did you use
> > to compute the magnitude? abs()? sqrt of product of complex number
> > with its conjugate?
> 
> I used the "0.0 + 0.0i" there to emphasize that I built PETSc with
> complex scalars, so all the scalars are stored as complex even if the
> imaginary components are always zero. I have been using std::abs()
> from the C++ standard library complex template class
> (std::complex<double>). I believe the values of type 'double' returned
> by std::abs() are implicitly cast to the real part of
> std::complex<double> or PetscScalar, though I usually try to cast
> explicitly.
> 
> > This seems like an application where the variables are of type
> > complex but the constraints need to be set for functions of the
> > variables which are actually real numbers. I am not sure whether
> > this can be supported. Maybe Barry or others can shed more
> > light on this.
> 
> My understanding is that comparison operators (such as < and >) are
> not defined for complex numbers, so would it be correct to say that
> anyone using complex scalars within PETSc will encounter this problem
> with SNESVI? In other words, SNESVI inequality constraints can only be
> applied to real variables or real-valued functions of complex
> variables?
> 
> >>
> >> The residuals are set as follows, and the Jacobian entries are set
> >> correspondingly:
> >> f[V_block] = {function of adjacent V_blocks, V_electrode}
> >> f[V_electrode] = 0.0 + 0.0i
> >> f[magnitude(V_electrode)] = 0.0 + 0.0i
> >> f[magnitude(I_electrode)] = 0.0 + 0.0i
> >
> > What equations do you use for the partial derivatives of the
> > magnitude functions w.r.t. the actual complex variables (V_block and
> > V_electrode) in the jacobian evaluation?
> 
> Since I set all the residuals for both magnitude variables to zero,
> and each entry in the Jacobian matrix is a partial derivative of a
> residual, I believe the Jacobian entries for both magnitudes would all
> be zero as well. However, your point made me realize I have no
> Jacobian entry for the partial derivative of f[V_block] w.r.t.
> V_electrode. Is that needed? How would one enter it with
> MatSetValuesStencil? My working understanding was that each degree of
> freedom conceptually has its own residual and Jacobian, even though
> all residuals are contained in one Vec and all Jacobians contained in
> one Mat.
> 
> >>
> >> The first two potentials, V_block and V_electrode, are independent,
> >> while magnitude(V_electrode) is a function of V_electrode only and
> >> the
> >> electrode current magnitude(I_electrode) is a function of V_block
> >> and
> >> V_electrode. Note that V_electrode acts as a source term in the
> >> finite
> >> difference formulation for V_block. While all of these variables
> >> are
> >> complex, the latter two always have zero imaginary parts. I suppose
> >> I
> >> am assuming that PETSc knows to compare only the real parts if the
> >> imaginary parts are zero. (Is this a bad assumption?)
> >>
> >> When solving with PETSc (using -snes_vi_type rs)
> >
> > Did you also use -snes_type vi?
> 
> I have been using SNESSetType(snes, SNESVI) in my code, and the use of
> -snes_type vi in the command line does not seem to change any
> behaviour.
> 
> >
> > V_electrode stays at
> >> its maximum, as set in the initial guess, and the V_block
> >> distribution
> >> falls properly from that. However, measurements of
> >> magnitude(I_electrode) are way above maximum and do not move
> >> downwards
> >> with successive SNES iterations. I can also set the constraint on
> >> magnitude(V_electrode) to less than maximum and it does not affect
> >> the
> >> value of V_electrode. How can I tell PETSc to change V_electrode
> >> when
> >> the magnitude(V_electrode) or magnitude(I_electrode) constraints
> >> are
> >> not met?
> >>
> >  Send the code for your application with instructions on how to run
> >  it and we'll try to figure out what needs to be
> >  done.
> 
> I really appreciate this generous offer of your time, and I will
> certainly take you up on it if we cannot resolve this otherwise.
> Unfortunately, I will be travelling for the next few weeks and it may
> take me some time to isolate the PETSc-dependent parts of my
> application. Please bear with me as I will try to keep in touch while
> I am away.
> 
> As I mentioned above, I think the issue may lie in the conceptual
> incompatibility of complex variables and inequality constraints, but I
> would appreciate your thoughts on that.
> 
> Thanks again,
> 
> Jon
> 
> >
> >
> > Shri
> >
> >
> >> Thanks again for your help.
> >>
> >> Jon
> >>
> >>
> >> On 2011-07-25, at 8:58 PM, Barry Smith wrote:
> >>
> >>>
> >>> On Jul 25, 2011, at 5:50 PM, Jonathan Backs wrote:
> >>>
> >>>> Hi Shri,
> >>>>
> >>>> Thanks for your message and all the helpful tips. If the
> >>>> TSVISetVariableBounds() functions are available now, I would like
> >>>> to try them as well. Is the interface essentially the same as
> >>>> with
> >>>> SNESVISetVariableBounds()? I will get back to you and Barry when
> >>>> I
> >>>> have had a chance to modify my application.
> >>>>
> >>>> For my problem, I believe it makes sense to have bounds on one of
> >>>> the variables as well as one function of the variables. The two
> >>>> relevant degrees of freedom are the block voltage (one for each
> >>>> finite difference block) and the electrode voltage (one for each
> >>>> electrode, which may be present in multiple blocks). The
> >>>> electrode
> >>>> voltage should keep a constant phase while the magnitude is
> >>>> constrained between zero and some maximum. The block voltages
> >>>> near
> >>>> the electrodes depend on the electrode voltages as well as the
> >>>> neighbouring block voltages. The electrode current for a given
> >>>> electrode is a function of its electrode voltage and several
> >>>> nearby
> >>>> block voltages, and should be constrained in magnitude between
> >>>> zero
> >>>> and some maximum (and the phase unconstrained). Would I need to
> >>>> use
> >>>> SNESVISetComputeVariableBounds since the electrode current is a
> >>>> function of the other variables? Would any other provisions need
> >>>> to
> >>>> be made for the block voltages, since they depend on the
> >>>> electrode
> >>>> voltages?
> >>>>
> >>>
> >>>  You cannot directly make a bound on some function of other
> >>>  variables. Instead you need to introduce another variable that is
> >>>  defined to be equal to that function and put the bound on that
> >>>  new
> >>>  variable.
> >>>
> >>>  Barry
> >>>
> >>>> Thank you again,
> >>>>
> >>>> Jon
> >>>>
> >>>> On 2011-07-25, at 11:58 AM, Shri wrote:
> >>>>
> >>>>>
> >>>>> ----- Original Message -----
> >>>>>> On Jul 22, 2011, at 4:16 PM, Jonathan Backs wrote:
> >>>>>>
> >>>>>>> Barry,
> >>>>>>>
> >>>>>>> Thank you so much for your response. Lucky, indeed! I look
> >>>>>>> forward
> >>>>>>> to trying out these new features.
> >>>>>>>
> >>>>>>> I neglected to mention in my original post that my electrical
> >>>>>>> problem is part of a DAE, which includes a time-dependent
> >>>>>>> heating
> >>>>>>> problem. Can SNESVI constraints be used in conjunction with
> >>>>>>> TSSetIFunction() and TSSetIJacobian() as well? (Of course, I
> >>>>>>> only
> >>>>>>> need the constraints for the time-independent electrical
> >>>>>>> portion.)
> >>>>>>
> >>>>>> We have not yet put that in but Shri is starting to look at
> >>>>>> that
> >>>>>> just
> >>>>>> now. Basically we would have a TSVISetVariableBounds() and
> >>>>>> handle
> >>>>>> everything from there. I suggest you start with a simple time
> >>>>>> electrical portion with constraints to explore and we'll go
> >>>>>> from
> >>>>>> there. Shri is actually a electrical networks guy and so can
> >>>>>> speak
> >>>>>> your language.
> >>>>>
> >>>>>
> >>>>> I've added TSVISetVariableBounds() for setting the bounds on the
> >>>>> variables using the TS object directly.
> >>>>> A few things that i want to mention here about using the
> >>>>> variational inequality nonlinear solver (SNESVI).
> >>>>> i) Use the runtime option -snes_type vi or explicitly set
> >>>>> SNESSetType(snes,SNESVI).
> >>>>> ii) There are two tested algorithms currently available, (a)
> >>>>> semismooth (-snes_vi_type ss) and (b) active set or reduced
> >>>>> space
> >>>>> (-snes_vi_type rs).
> >>>>> iii) Take a look at example,ex61.c, in
> >>>>> src/snes/examples/tutorials
> >>>>> which is a time-stepping nonlinear problem with constraints on
> >>>>> the
> >>>>> variables. This example does not use the TS object directly but
> >>>>> rather a time-loop is explicitly written. Another example,ex8.c,
> >>>>> in src/snes/examples/tests/ is a minimum surface area problem
> >>>>> which uses SNESVI.
> >>>>>
> >>>>> By the way, does your problem have bounds on the variables or
> >>>>> bounds on some function of the variables?
> >>>>>
> >>>>> Shri
> >>>>>
> >>>>>
> >>>>>
> >>>>>>
> >>>>>> Barry
> >>>>>>
> >>>>>>
> >>>>>>
> >>>>>>>
> >>>>>>> Thank you again,
> >>>>>>>
> >>>>>>> Jon
> >>>>>>>
> >>>>>>> On 2011-07-22, at 2:46 PM, Barry Smith wrote:
> >>>>>>>
> >>>>>>>>
> >>>>>>>> Jon,
> >>>>>>>>
> >>>>>>>> You may be in luck. In PETSc-dev
> >>>>>>>> http://www.mcs.anl.gov/petsc/petsc-as/developers/index.html
> >>>>>>>> we
> >>>>>>>> have now implemented variational inequality non-linear
> >>>>>>>> solvers
> >>>>>>>> with box constraints.
> >>>>>>>>
> >>>>>>>> That is one has a set of variables u and algebraic equations
> >>>>>>>> F(u) = 0 (say of size n each) but in addition one has
> >>>>>>>> constraints lu <= u <= uu where some or all of lu may be
> >>>>>>>> negative infinity (no constraints) and some or all of uu may
> >>>>>>>> be
> >>>>>>>> infinity (no constraints). There is also a constraint on the
> >>>>>>>> sign of F() for those equations associated with active
> >>>>>>>> constraints. If your constraints are not in this form
> >>>>>>>> sometimes
> >>>>>>>> you can introduce new additional variables to get it in this
> >>>>>>>> form. Read up a little on variational inequalities on the
> >>>>>>>> web.
> >>>>>>>>
> >>>>>>>> To use this you provide the usual SNES function and Jacobian
> >>>>>>>> but
> >>>>>>>> you also provide SNESVISetVariableBounds() there is a manual
> >>>>>>>> page for this function and for for SNESVI at
> >>>>>>>> http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-dev/docs/index.html
> >>>>>>>> (the link is broken right now but Satish is fixing it). Plus
> >>>>>>>> several examples in src/snes/examples/tutorials (in
> >>>>>>>> petsc-dev).
> >>>>>>>>
> >>>>>>>> This is all new code so we would be interested in your
> >>>>>>>> feedback.
> >>>>>>>>
> >>>>>>>>
> >>>>>>>> Barry
> >>>>>>>>
> >>>>>>>>
> >>>>>>>>
> >>>>>>>>
> >>>>>>>>
> >>>>>>>> On Jul 22, 2011, at 3:33 PM, Jonathan Backs wrote:
> >>>>>>>>
> >>>>>>>>> Hi,
> >>>>>>>>>
> >>>>>>>>> I am trying to add a constraint feature to my first PETSc
> >>>>>>>>> application, which uses the finite difference method to
> >>>>>>>>> calculate
> >>>>>>>>> the potential distribution produced by a collection of
> >>>>>>>>> electrodes
> >>>>>>>>> in a resistive medium. I would like to make this simulation
> >>>>>>>>> more
> >>>>>>>>> realistic by imposing a maximum electric current and a
> >>>>>>>>> maximum
> >>>>>>>>> potential difference that can be supplied to each electrode
> >>>>>>>>> by
> >>>>>>>>> the
> >>>>>>>>> power supply. If the medium between the electrodes is very
> >>>>>>>>> conductive, the current maximum would be exceeded by the
> >>>>>>>>> maximum
> >>>>>>>>> potential difference, so the potential difference should be
> >>>>>>>>> decreased from maximum until it produces the maximum
> >>>>>>>>> current.
> >>>>>>>>> On
> >>>>>>>>> the other hand, the potential difference between the
> >>>>>>>>> electrodes
> >>>>>>>>> should remain at maximum as long as the current remains
> >>>>>>>>> below
> >>>>>>>>> maximum (say, for a less conductive medium).
> >>>>>>>>>
> >>>>>>>>> I added an extra degree of freedom (the electrode voltages)
> >>>>>>>>> to
> >>>>>>>>> my
> >>>>>>>>> DMDA, and I developed a set of conditional expressions that
> >>>>>>>>> describe the above constraints, but one problem is that the
> >>>>>>>>> logic
> >>>>>>>>> relies on if-then-else decisions that are made when forming
> >>>>>>>>> the
> >>>>>>>>> function/residual and the Jacobian. Once these decisions are
> >>>>>>>>> made,
> >>>>>>>>> of course, the conditions are not checked again until the
> >>>>>>>>> next
> >>>>>>>>> function or Jacobian evaluation. The non-linear solver then
> >>>>>>>>> tends
> >>>>>>>>> to oscillate between extreme solutions to the opposing
> >>>>>>>>> conditions
> >>>>>>>>> with each iteration, and never converges towards a
> >>>>>>>>> reasonable
> >>>>>>>>> solution.
> >>>>>>>>>
> >>>>>>>>> Is there a better strategy for solving such problems? Does
> >>>>>>>>> PETSc
> >>>>>>>>> offer mechanisms to aid in their solution? I would very much
> >>>>>>>>> appreciate any hints.
> >>>>>>>>>
> >>>>>>>>> Thank you for your time,
> >>>>>>>>>
> >>>>>>>>> Jon
> >>>>>>>>
> >>>>>>>
> >>>>>
> >>>>
> >>>
> >