[petsc-users] Conditional Constraints

[Jonathan Backs]

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