[petsc-users] Conditional Constraints

[Jonathan Backs]

Shri,

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

Would f(Ve_ijk) = Imag - Imag_max = 0 simply be equivalent to saying, solve the real system of equations Imag - Imag_max = 0 and 0 = 0 (the imaginary part)? Regardless, I am not sure what a correct choice of f(Ve_ijk) would be.


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


Yes, that is essentially what I am trying to find. Ve determines Vb, and the difference between the two determines Imag. Vmag (magnitude of Ve) and Imag have the restrictions 0 < Vmag < Vmag_max and 0 < Imag < Imag_max. I want to find either Ve for Imag = Imag_max or Imag for Vmag = Vmag_max. In my test scenario, Vmag_max happens to be quite high so Imag_max will always be reached first, but in general either situation is possible.

Thank you again,

Jon


On 2011-09-08, at 8:34 PM, Shri wrote:

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

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