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>
