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?
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
>
> 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/9cc37bc1/attachment-0001.htm>
