[petsc-users] Conditional Constraints

Thu Sep 8 21:34:55 CDT 2011

Jon, 

> I do not want Ve to be constant, however; I would like it to either decrease in response to Imag > Imag_max or increase in response to Imag < Imag_max, providing > the condition Vmag <= Vmag_max remains true. 
> Should I then have something like f(Ve_ijk) = Imag - Imag_max? 

Note here that (a) f(Ve_ijk) is a complex function while Imag and Imag_max are real. So the above equation does not hold. 
(b) In your previous email, you had bounds on Imag, 0 < Imag < Imag_max. Now you want Ve_ijk to increase if Imag > Imag_max and decrease if Imag < Imag_max. Are you trying to find the value of Ve at which the current is maximum, i.e. Imag = Imag_max, by doing this? 

----- Original Message -----

Shri, 

Thanks very much for your response. In fact, -snes_converged_reason revealed that the solution was not converging, as you suggested: "Nonlinear solve did not converge due to DIVERGED_LINE_SEARCH." I had not realized this earlier since the voltage distribution in the solution appeared to be correct, even though the SNES_VI constraints were not satisfied. 

Given that and your explanation, I am sure something is wrong with my choices of f. I do not want Ve to be constant, however; I would like it to either decrease in response to Imag > Imag_max or increase in response to Imag < Imag_max, providing the condition Vmag <= Vmag_max remains true. Should I then have something like f(Ve_ijk) = Imag - Imag_max? 

Note that you have 

One more thing from your response: 

>If a block (with coordinates i, j, k) is not occupied by an electrode , C_e is set to zero in the expression for Vb_ijk and all function and Jacobian entries for Ve, >Vmag, and Imag are either not set or set to zero. 

Does this mean f(Ve_ijk) = f(Vmag_ijk) = f(Imag_ijk) = 0? If so, then Ve_ijk, Vmag_ijk, and Imag_ijk would retain their initial values for all time steps. 

In my present scheme, f(Ve_ijk) = f(Vmag_ijk) = 0 and f(Imag_ijk) is non-zero for electrode blocks; all three are zero for non-electrode blocks. Depending on the choice of initial conditions, I can make f(Vmag_ijk) non-zero as well (until Ve changes to match Ve_max). 

Thanks again, 

Jon 

On 2011-09-08, at 3:42 PM, Shri wrote: 

Jon, 

> for Ve_ijk, f = 0.0 + 0.0*i, and all corresponding Jacobian entries are zero as well. Ve_ijk is set to the maximum voltage for each electrode when the initial solution >is - set. Since the magnitude of Ve is constrained rather than Ve itself, I believe this is correct, but could this be why Ve does not change as it should when the >problem is - being solved? 

This f(Ve_ijk) seems very odd to me since normally you would always have f(some/all variables) = 0 as the function. If f=0.0+0.0*i for all steps, then you don't need a seperate variable Ve_ijk, you can treat it as a constant since Ve_ijk will not change its value during the solutions. 

Moroever, the Jacobian rows for df(Ve_ijk)/d(all variables) would be all zeros which implies a singular matrix. I am not sure how the solution converged for this case. Were you using -pc_type lu with some shift or some external direct solver such as superlu? What do you get for -snes_monitor -snes_converged_reason? 

> Ve_ijk is set to the maximum voltage for each electrode when the initial solution is set. Since the magnitude of Ve is constrained rather than Ve itself. 
> for Vmag_ijk, f = abs(Ve_ijk) - abs(Ve_ijk_max) 

Since Ve_ijk = Ve_ijk_max, f for Vmag_ijk = 0.0 and hence Vmag_ijk also will also retain its initial value throughout. 

>If a block (with coordinates i, j, k) is not occupied by an electrode , C_e is set to zero in the expression for Vb_ijk and all function and Jacobian entries for Ve, >Vmag, and Imag are either not set or set to zero. 
Does this mean f(Ve_ijk) = f(Vmag_ijk) = f(Imag_ijk) = 0? If so, then Ve_ijk, Vmag_ijk, and Imag_ijk would retain their initial values for all time steps. 

> When I run my program with -snes_vi_type rs, the solution converges quickly but Vmag voltages remain the same as Ve_max for every electrode. As a result, Imag > currents are far larger than abs(iE_max). 
Seems to me like the problem is not converging, or giving the incorrect solution, since the solution of a variational inequality problem has to satisfy Imag < iE_max, which in your case is not happening. 

Thanks for the detailed explanation. Would you mind if I (or you) forward this to petsc-users list so that others could share their thoughts too. 

Shri 

----- Original Message -----

Hi Shri, 

Thanks again for your patience with my being away this past month. I have spent some more time with this problem and I have added Jacobian entries for 'cross terms' such as d(f(Vb)) / d(Ve). However, this did not appear to change the solution arrived at by PETSc. 

I will give you some examples of the f functions and corresponding Jacobian entries as you suggested, with the hope you will be able to see where I may be going wrong. I will use Vb = V_block, Ve = V_electrode, Vmag = magnitude(V_electrode), Imag = magnitude(I_electrode) to name the four degrees of freedom. 

In my 3D finite differencing grid, there are two types of blocks: blocks occupied by electrodes, and blocks occupied by a resistive medium. The block type of the current block and its neighbouring blocks determines the form of f. This is complicated by the fact that individual electrodes generally occupy three vertically-contiguous blocks, and I want to limit the electric current flowing from the whole electrode rather than the current through each block. Note that Ve, Vmag, and Imag are thus computed or recorded only for the top block of each electrode. 

If a block (with coordinates i, j, k in the x, y, and z directions, respectively) is occupied by an electrode , 
- for Vb_ijk, f = C_posX * Vb_posX + C_posY * Vb_posY + C_posZ * Vb_posZ 
- (C_posX + C_posY + C_posZ + C_negX + C_negY + C_negZ - C_e) * Vb_ijk 
+ C_negX * Vb_negX + C_negY * Vb_negY + C_negZ * Vb_negZ 
- C_e * Ve, 
where all C's are calculated constants that do not depend on the other degrees of freedom, but do depend on the type of the neighbouring block in their assigned directions. Then, 
the Jacobian entries for Vb_ijk are 
df/dVb_negZ = C_negZ (col.i = i, col.j = j, col.k = k - 1) 
df/dVb_negY = C_negY (col.i = i, col.j = j - 1, col.k = k) 
df/dVb_negX = C_negX (col.i = i - 1, col.j = j, col.k = k) 
df/dVb_ijk = - (C_posX + C_posY + C_posZ + C_negX + C_negY + C_negZ - C_e) (col.i = i, col.j = j, col.k = k) 
df/dVb_posX = C_posX (col.i = i + 1, col.j = j, col.k = k) 
df/dVb_posY = C_posY (col.i = i, col.j = j + 1, col.k = k) 
df/dVb_posZ = C_posZ (col.i = i, col.j = j, col.k = k + 1) 
For each of these, row.i = i, row.j = j, row.k = k, and row.c = DOF_VB (the DOF index for Vb). Also, col.c = DOF_VB. 
The only cross term for the Vb_ijk function is 
df/dVe_ijk = -C_e (row.i = i, row.j = j, row.k = k, row.c = DOF_VB, col.i = i, col.j = j, col.k = k col.c = DOF_VE) 

- for Ve_ijk, f = 0.0 + 0.0*i, and all corresponding Jacobian entries are zero as well. Ve_ijk is set to the maximum voltage for each electrode when the initial solution is set. Since the magnitude of Ve is constrained rather than Ve itself, I believe this is correct, but could this be why Ve does not change as it should when the problem is being solved? 

- for Vmag_ijk, f = abs(Ve_ijk) - abs(Ve_ijk_max). The idea was to provide some 'pressure' towards Ve = Ve_max. The corresponding Jacobian entries are 
df/dVmag_ijk = 1.0 + 0.0*i (row.i = i, row.j = j, row.k = k, row.c = DOF_VMAG, col.i = i, col.j = j, col.k = k, col.c = DOF_VMAG) 
df/dVe_ijk = 0.5 * conj(Ve_ijk) / abs(Ve_ijk) (row.i = i, row.j = j, row.k = k, row.c = DOF_VMAG, col.i = i, col.j = j, col.k = k, col.c = DOF_VE) 
(However, since d(abs(Ve_ijk))/dVe_ijk does not have satisfied Cauchy-Reimann equations for the entire complex plane, this derivative does not technically exist.) 

- for Imag_ijk, f = abs(iE) - iE_max. Again, the idea was to provide some 'pressure' towards iE = iE_max. iE is calculated by adding iE_block currents for each block occupied by an electrode. iE_block for block i,j,k is found from 
iE_block_ijk = aE_ijk * (Ve - Vb_ijk) 
where aE_ijk is a constant that does not depend on any degrees of freedom. There is one Ve value per electrode. The Jacobian entries are then 
df/dImag_ijk = 1.0 + 0.0*i (row.i = i, row.j = j, row.k = k, row.c = DOF_IMAG, col.i = i, col.j = j, col.k = k, col.c = DOF_IMAG) 
df/dVe = 0.5 * aE_sum * conj(iE) / abs(iE) (row.i = i, row.j = j, row.k = k, row.c = DOF_IMAG, col.i = i, col.j = j, col.k = k, col.c = DOF_VE) 
df/dVb_ijk = - 0.5 * aE_ijk * conj(iE) / abs(iE) (row.i = i, row.j = j, row.k = k, row.c = DOF_IMAG, col.i = i, col.j = j, col.k = k, col.c = DOF_VB) 
df/dVb_ij(k-1) = - 0.5 * aE_ij(k-1) * conj(iE) / abs(iE) (row.i = i, row.j = j, row.k = k, row.c = DOF_IMAG, col.i = i, col.j = j, col.k = k-1, col.c = DOF_VB) 
df/dVb_ij(k-2) = - 0.5 * aE_ij(k-2) * conj(iE) / abs(iE) (row.i = i, row.j = j, row.k = k, row.c = DOF_IMAG, col.i = i, col.j = j, col.k = k-2, col.c = DOF_VB) 

If a block (with coordinates i, j, k) is not occupied by an electrode , C_e is set to zero in the expression for Vb_ijk and all function and Jacobian entries for Ve, Vmag, and Imag are either not set or set to zero. 

My SNES_VI constraints are set as follows for every block: 
xl_ijk_DOF_V = -SNES_VI_INF 
xu_ijk_DOF_V = SNES_VI_INF 

xl_ijk_DOF_VE = -SNES_VI_INF 
xu_ijk_DOF_VE = SNES_VI_INF 

xl_ijk_DOF_VMAG = 0.0 
xu_ijk_DOF_VMAG = abs(Ve_max) 

xl_ijk_DOF_IMAG = 0.0 
xu_ijk_DOF_IMAG = abs(iE_max) 

When I run my program with -snes_vi_type rs, the solution converges quickly but Vmag voltages remain the same as Ve_max for every electrode. As a result, Imag currents are far larger than abs(iE_max). Are you able to see from here what the problem might be? Again, my hunches are that the complex derivatives do not work in my Jacobian because they do not exist, or that my choices of f for Ve, Vmag, and Imag are not correct. 

Please let me know if you would like clarification on any portion of my description., and thank you again for all your help with this. 

Jon 

On 2011-07-29, at 12:13 PM, Shri wrote: 

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 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20110908/0af5c3db/attachment-0001.htm>