[Getdp] Multiple right hand side for the same stiffness matrix

[Benjamin Vial]

Thanks Christophe, I will have a look at this.


Yes I remember that at some point I saw some codes about optimization

in the sources of getp/gmsh/onelab (can't remember)


Thanks again.


Kind regards,


Ben


---------------------------------------------------------------------------------------------------
Benjamin Vial
Postdoctoral Research Assistant
Antennas and Electromagnetics research group
School of Electronics Engineering and Computer Science
Queen Mary, University of London
http://benjaminvial.org/
http://antennas.eecs.qmul.ac.uk/
http://www.quest-spatial-transformation.org/team/researchers/dr-benjamin-vial/

________________________________
From: Christophe Geuzaine <cgeuzaine at ulg.ac.be>
Sent: 09 September 2016 09:33:03
To: Benjamin Vial
Cc: getdp at onelab.info; Erin Kuci
Subject: Re: [Getdp] Multiple right hand side for the same stiffness matrix


> On 08 Sep 2016, at 20:07, Benjamin Vial <b.vial at qmul.ac.uk> wrote:
>
> Dear all,
>
> First I would like to thank all getdp and gmsh authors/developers for such amazing pieces of software
> and to make them available open source!
>
> I am currently working on topology optimization problems in electromagnetism.
> I'll try to keep it simple: I want to maximize square the norm of the electric field in some sub-domain Omega_target,
> by changing the permittivity value in a design sub-domain Omega_design.
> I defined a formulation "helmoltz_scalar" for a 2D time-harmonic problem (xy plane), TE polarization
> for the z component of the E-field, which I call u.
> In matrix form :
> A*u = b   (1)
> where b is a known source term.
>
> Now gradient-based optimization algorithm needs information about the sensitivity of the
> objective function g with respect to the design variable epsilon. The most efficient way to obtain it
> is by solving an adjoint problem, namely
> A*uadj = dg/du     (2)
> In my case g=|u|^2 so the second member of the adjoint equation (2) is dg/du = 2*Re[u]
>
> Since the matrix A is the same for both problems, we do not need to assemble it twice
> Is there a way to save time and generate only the right hand side for the adjoint equation?
> It looks like the function GenerateRHSGroup might be useful (but it is not documented yet), together with SolveAgain.
>

Indeed, these are the two functions you will need: SolveAgain will reuse the factorization (or in general the preconditioner) if it's available; and GenerateRHSGroup[A, region] will recompute the right-hand-side, and only evaluate the equation terms defined on "region".

I'm CC:ing Erin Kuci, who is developing a Python optimization framework for ONELAB, for both shape and topology optimization.

Christophe

> This somehow similar to this old thread but maybe there is a more straightforward way to do it?
> http://onelab.info/pipermail/getdp/2011/001383.html
>
>
> Here are the Formulation and Resolution parts of my .pro file:
>
> #############################################################################
> Formulation{
>     {Name helmoltz_scalar; Type FemEquation;
>         Quantity {
>         { Name u; Type Local; NameOfSpace Hgradu;}
>         }
>         Equation {
>         Galerkin { [k0^2*CompZZ[epsilonr[]]*Dof{u} , {u}];
>         In Omega; Jacobian JVol; Integration Int_1;  }
>         Galerkin { [-1/TensorDiag[CompYY[mur[]],CompXX[mur[]],CompXX[mur[]]]*Dof{Grad u} , {Grad u}];
>         In Omega; Jacobian JVol; Integration Int_1; }
>         Galerkin { [source[] , {u}];
>         In Omega; Jacobian JVol; Integration Int_1;  }
>         }
>     }
>     {Name adjoint_eq; Type FemEquation;
>         Quantity {
>         { Name u; Type Local; NameOfSpace Hgradu;}
>         { Name uadj; Type Local; NameOfSpace Hgraduadj;}
>         }
>         Equation {
>         Galerkin { [  k0^2*CompZZ[epsilonr[]]*Dof{uadj} , {uadj}];
>         In Omega; Jacobian JVol; Integration Int_1;  }
>         Galerkin { [-1/TensorDiag[CompYY[mur[]],CompXX[mur[]],CompXX[mur[]]]*Dof{Grad uadj} , {Grad uadj}];
>         In Omega; Jacobian JVol; Integration Int_1; }
>         Galerkin { [2*Re[{u}], {uadj}];  //
>         In Omega_target; Jacobian JVol; Integration Int_1;  }
>         }
>     }
>
> }
>
>
> Resolution {
>     { Name all
>         System {
>         { Name S; NameOfFormulation helmoltz_scalar; Type ComplexValue; Frequency Freq;}
>         { Name Sadj; NameOfFormulation adjoint_eq; Type ComplexValue; Frequency Freq;}
>         }
>         Operation {
>         Generate[S] ;Solve[S] ;SaveSolution[S] ;
>         Generate[Sadj] ;Solve[Sadj] ; SaveSolution[Sadj] ;
>         }
>
>     }
> }
>
> #############################################################################
>
>
>
> Thanks in advance!
>
> All best
>
> Benjamin
>
>
> ---------------------------------------------------------------------------------------------------
> Benjamin Vial
> Postdoctoral Research Assistant
> Antennas and Electromagnetics research group
> School of Electronics Engineering and Computer Science
> Queen Mary, University of London
> http://benjaminvial.org/
> http://antennas.eecs.qmul.ac.uk/
> http://www.quest-spatial-transformation.org/team/researchers/dr-benjamin-vial/
>
> _______________________________________________
> getdp mailing list
> getdp at onelab.info
> http://onelab.info/mailman/listinfo/getdp

--
Prof. Christophe Geuzaine
University of Liege, Electrical Engineering and Computer Science
http://www.montefiore.ulg.ac.be/~geuzaine

Free software: http://gmsh.info | http://getdp.info | http://onelab.info

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://onelab.info/pipermail/getdp/attachments/20160914/d9d529b6/attachment.html>