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

[Christophe Geuzaine]

> 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