Hi Shri,
I have continued to have little success with my SNESVI formulation of the current constraint. However, working outside of PETSc, I believe I have found an explicit formulation that works well:
1) Set Vmag = Vmag_max and use PETSc to solve for the electric potential distribution (V_ijk is the only degree of freedom).
2) For each electrode, let Vmag_new = Imag_max * Vmag_old / Imag_old.
3) Repeat until Imag is sufficiently close to Imag_max, or until Vmag = Vmag_max with Imag < Imag_max.
This process converges in about 5 iterations, which is great, but I would like to incorporate the thermal problem again (including TS) now that this electrical problem is solved. How could I implement this constraint inside the PETSc solution loop? Is there any advantage to using SNESVI when my update step is explicit? Or is there a feature I could use to add a Vmag_new computation after each SNES step, and have this taken into account for the SNES convergence criteria?
Thanks again for your continued patience and interest in helping me with this problem.
Jon
On 2011-09-15, at 12:07 PM, Shri wrote:
> Jon,
>
> > Additionally, the potential distribution does not look correct as it did before (even without SNES reporting convergence)---instead, the magnitude of the distribution > looks like the real part of the correct distribution, if that makes sense.
> Hmm, in the SNESVI formulation using complex numbers, we have the check real(xlower) <= real(x) <= real(xupper). i am not sure whether this is the cause for such a behavior. I've had conversation with the other PETSc developers to use absolute values for comparison instead of real values, i.e. abs(xl) <= abs(x) <= abs(xupper). By doing this, you can essentially have a function for Ve and its bounds being 0 and Vmag_max. However, the other developers think that absolute values being nonlinear in nature might produce some odd behavior in the solution process, so i've havent' support for it yet.
>
> >This asymmetry may be due to the fact that I now only record the Vmag and Imag variables for the top block of each electrode, whereas before I recorded Ve for all >three blocks.
> Yes, do you want the magnitude and angle for all the three blocks to be same? If so, then probably you'll have to add another constraint.
>
> >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 Vmag >and Imag are either not set or set to zero.
> Do you already know the Vmag and Imag values are these blocks? Are they zero/non-zero? I suggest the following for these blocks
>
> xl_ijk_DOF_VMAG = Vmag0_ijk
> xu_ijk_DOF_VMAG = Vmag0_ijk
> xl_ijk_DOF_IMAG = Imag0_ijk
> xu_ijk_DOF_IMAG = Imag0_ijk
>
> where Vmag0_ijk and Imag0_ijk is the initial voltage and current magnitudes at these blocks. Basically, the idea here is to keep these magnitudes constant and not have the functions for these variables be a part of the solution process.
>
> Shri
>
>
>
>
>
> Hi Shri,
>
> Thanks for following up with me. I have made some progress, although my problem is not yet solved.
>
> Inspired by your recent suggestions, I realized Ve was redundant as a degree of freedom and did not need to be computed separately from Vmag. I do not want the phase of Ve to be modified by the solution process anyways, so whenever I need the complex value of Ve, I can simply compose it from the Vmag and the constant phase I have recorded for a given electrode. My new function and Jacobian formulations are below if you are interested in seeing them.
>
> After I made this change, I began to see the desired behaviour where Vmag decreases when Imag is too high. However, it seems to converge on some magnitude that is maybe 200% too high, and the SNES reports DIVERGED_LINE_SEARCH. Additionally, the potential distribution does not look correct as it did before (even without SNES reporting convergence)---instead, the magnitude of the distribution looks like the real part of the correct distribution, if that makes sense.
>
> A further problem relates to the fact that each electrode is three blocks tall. Before I eliminated Ve, each of the three blocks had a very similar voltage and phase. Now, the top block has a different magnitude and phase than the two beneath it. This asymmetry may be due to the fact that I now only record the Vmag and Imag variables for the top block of each electrode, whereas before I recorded Ve for all three blocks. As you can find below, however, I have included Jacobian entries for the lower two block voltages Vb of each electrode, so I cannot see why this strategy would not work.
>
> Given your experience with PETSc and Newton's method, do you have any ideas about where these problems might originate?
>
> Thank you again,
>
> Jon
>
>
> 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 * polar(Vmag, arg(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/dVmag = polar(-C_e, arg(Ve)) (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_VMAG)
>
> - 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)
>
> - 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/dVmag = aE_sum * (real(iE) * cos(arg(Ve)) + imag(iE) * sin(arg(Ve)))/ 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_VMAG)
> 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)
> (These last three are the only derivatives left that do not technically exist.)
>
> 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 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_VMAG = 0.0
> xu_ijk_DOF_VMAG = abs(Ve_max)
> xl_ijk_DOF_IMAG = 0.0
> xu_ijk_DOF_IMAG = abs(iE_max)
>
>
> On 2011-09-13, at 12:12 PM, Shri wrote:
>
> Jon,
> Any progress with your work? Send us an email once you've sorted out your formulation and we would be glad to help you out with its PETSc implementation.
>
> Shri
>
> 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
>
>
>
>
>
>
