> 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)?
Yes.
> Regardless, I am not sure what a correct choice of f(Ve_ijk) would be
Maybe, having the equation for current block
iE_block_ijk = aE_ijk * (Ve - Vb_ijk) for f(Ve_ijk)?
where iE_block_ijk = Imag_max + j*0 ?
Perhaps, you can have a different and better expression for f(Ve_ijk) based on
some related physics.
----- Original Message -----
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?
----- 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL:
<http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20110909/8c9eaf8a/attachment-0001.htm>
