[petsc-users] Conditional Constraints
Jonathan Backs
jbacks at mcmillan-mcgee.com
Thu Sep 8 18:07:57 CDT 2011
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>
More information about the petsc-users
mailing list