[Nek5000-users] Clarification on Adjoint solver

Wed Feb 10 11:58:34 CST 2016

Adam and Oana,
Thanks again for your responses.
Let me summarize what I understood.
1) We may have either evolving (unsteady) or non-evolving (steady) base 
flow. As for the latter, we may use restart option in rea file.
2) For non-evolving case (which is what I have in mind in the end) Nek 
"is able to" handle and to do so parameter #31 in rea file needs to be 
modified. Can you please clarify how. Suppose that there is velocity and 
a passive scalar or even temperature, so there are three adjoint fields 
(velocity, pressure and the scalar).
Also, is what Oana mentioned, "one has to do the change of the 
velocities himself" related to remark number 2 above?
3) There's also a confusion in userchk: First Adam mentioned that we 
need to include 'ADJOINT' plus if(ISTEP.eq.0) ADJOINT=.TRUE. for adjoint 
solver in order to i) changing the time and pressure sign and 2) 
modification of advection terms. But in the second example, while he 
still uses include 'ADJOINT' in userchk but then there is this statement 
IFADJ = .FALSE. Can you please clarify what's the difference here?

I'm really grateful that you take the time to answer these questions. 
Please let me know if we can take this conversation off the mailing list.

Brief summary of the problem we have in mind: we'd like to optimize the 
flow with respect to some parameters (either BCs or ICs) in terms of a 
defined cost function. So this can be regarded as either design or 
optimal control problem. Calculus of variations and Lagrange multipliers 
lead to adjoint equations which are very similar to NS equations. We 
already solved the problem using NEK and we'd like to use this again for 
adjoint equations. All the results will be in an outer loop for 
optimization (based on, say, steepest gradient). In general, the problem 
is unsteady, but we can start off with steady and non-evolving cases at 
early stages.

Best,
Saleh

On 2/8/2016 4:26 AM, nek5000-users at lists.mcs.anl.gov wrote:
> Hi Saleh
>
> If I understand you correctly, you ask about a base flow, i.e. the 
> velocity field you do your linearisation about (the one you use to 
> advect prturbation). If this is the case, than it is possible to do 
> with nek. Otherwise you have to wait for the code Oana mentioned. What 
> we do in our stability calculations, we take the solution of 
> non-linear NS (steady state reached with e.g. SFD) and next we load it 
> into nek as non-linear filed. One can do it by setting restart options 
> in ###.rea file. It should look like:
>    3.33333       3.33333     -0.833333      -1.16667 
> XFAC,YFAC,XZERO,YZERO
>  **MESH DATA** 1st line is X of corner 1,2,3,4. 2nd line is Y.
>        -1472  2          1472     NELT,NDIM,NELV
>            1  PRESOLVE/RESTART OPTIONS  *****
> bfext_cyl?.f00000
>            7          INITIAL CONDITIONS *****
> C Default
> C Default
>
>
> Where bfext_cyl0.f00000 is our restart file (we usually save this 
> velocity field in double precision). Next you have to set linear 
> solver with non-evolving base flow (nek is able to evolve both so be 
> careful). You freeze baseflow by setting negative value of parameter 
> 31 in ###.rea:
>    0.00000     p030 > 0 ==> properties set in uservp()
>   -1.00000     p031 NPERT: #perturbation modes
>    0.00000     p032 #BCs in re2 file, if > 0
>
> Nek is able of evolving number of perturbation fields set by 
> abs(param(31)) (in ###.rea) and lpert in SIZE. I have to mention that 
> nek by itself outputs the non-linear fields only, so to see 
> perturbation you should add to your userchk something like this
>
>
>       subroutine userchk
>
>       include 'SIZE'            ! NX1, NY1, NZ1, NELV, NIO
>       include 'TSTEP'           ! ISTEP, IOSTEP
>       include 'SOLN'            ! V[XYZ]P, TP
>       include 'INPUT'           ! IF3D, PARAM
>       include 'ADJOINT'         ! IFADJ
>
> !     direct or adjoint run
>       IFADJ = .FALSE.
>
> !     write perturbation field
>       if (mod(ISTEP,int(PARAM(15))).eq.0) then
>          call outpost2(VXP,VYP,VZP,PRP,TP,0,'prt')
>       endif
>
> I hope it is an answer to your question.
> regards
> adam
>
> On 2016-02-06 19:15, nek5000-users at lists.mcs.anl.gov wrote:
>> Hi Saleh,
>>
>> There is no easily accessible implementation for this, one has to do 
>> the change of the velocities himself. In fact we have ongoing work on 
>> this and once we are done I guess we can share it.
>> What kind of case are you looking at, long time horizons are quite 
>> tricky to handle. We can take this conversation off the mailing list 
>> if needed.
>>
>> Oana
>>
>>
>>
>> ________________________________________
>> From: nek5000-users-bounces at lists.mcs.anl.gov 
>> [nek5000-users-bounces at lists.mcs.anl.gov] on behalf of 
>> nek5000-users at lists.mcs.anl.gov [nek5000-users at lists.mcs.anl.gov]
>> Sent: Friday, February 05, 2016 2:15 PM
>> To: nek5000-users at lists.mcs.anl.gov
>> Subject: Re: [Nek5000-users] Clarification on Adjoint solver
>>
>> Adams,
>>
>> I really appreciate your detailed answer. I definitely clarified the Nek
>> functionality for adjoint solvers. One question remains is how to
>> include the velocity of the forward solver into adjoint equations. I
>> assume this is also the case when you do linear stability analysis with
>> Nek, since the velocity in convective terms is not completely part of
>> the solution. Mathematically, my question is:
>>
>> in $V. \nabla u_a$ term, as the advection, how can we substitute $V$ (as
>> the forward velocity which is already solved for in a regular NS
>> solution) in adjoint equations. suppose that the velocity, along with
>> other quantities are saved in afield out put such as fld.f0000X for X
>> time-step. So we have access to that velocity but how can we include it
>> in adjoint equations?
>>
>> Many thanks again,
>> Saleh
>>
>> On 2/2/2016 3:09 AM, nek5000-users at lists.mcs.anl.gov wrote:
>>> Hi
>>>
>>> The adjoint implementation in nek assumes change of the direction of
>>> time evolution. The goal is to use the standard nek solver with
>>> minimal modifications (only advection terms have to be changed). You
>>> get these equations from you standard one by changing the time and
>>> pressure sign (t -> -t; p-> -p). This of course changes your time
>>> integration limits, but it shouldn't be a problem. So when you  have
>>> your linear direct simulation, to turn it into dual one it is enough
>>> to include  in your userchk include files ADJOINT and TSTEP, and for
>>> ISTEP.eq.0 set IFADJOINT to true. Something like:
>>>        subroutine userchk
>>>
>>>        include 'SIZE'
>>>        include 'TSTEP'           ! ISTEP
>>>        include 'ADJOINT'            ! IFADJ
>>>
>>>        if(ISTEP.eq.0) ADJOINT=.TRUE.
>>>
>>> There is no need to modify ###.rea or ###.map files. The only problem
>>> are boundary conditions for open flows, as they would require
>>> non-homogeneous bc. In this case to avoid problems we usually use
>>> sufficiently large domain together with zero Dirichlet bc. The last
>>> issue is sufficient resolution for both direct and adjoin runs. I hope
>>> this gives you some idea about nek implementation.
>>> Regards
>>> Adam
>>>
>>> On 2016-02-01 22:17, nek5000-users at lists.mcs.anl.gov wrote:
>>>> Hi Neks,
>>>>
>>>> To my knowledge, Nek5000 is able to solve adjoint equations of the 
>>>> form:
>>>> $\partial u_a/\partial t + V.\nabla u_a - u_a.(\nabla V)^T + \nabla
>>>> \p_a + 1/Re \nabla^2 u_a=0$ and $\nabla.u_a=0$
>>>> where u_a and p_a are adjoint velocity and pressure. (Slightly
>>>> different formulation may be seen in literature).
>>>>
>>>> Let's assume the forward problem is already solved so that $V$, i.e.
>>>> velocity, (and possibly $p$ as the pressure) is (are) already know.
>>>> Thus, the initial and boundary conditions for adjoint NS are also
>>>> known. How can we now solve the dual/adjoint NS problem with Nek5000?
>>>> Is there any example on how to modify .rea files? Specailly, the
>>>> "convection term of $V.\nabla u_a - u_a.(\nabla V)^T$ needs to be
>>>> modified compared to forward NS but am not sure how.
>>>> My search of mailing list entails some modifications in "perturb.f"
>>>> subroutine, but a little bit more clarification on that would be
>>>> really appreciated.
>>>>
>>>> Yours,
>>>> Saleh
>>>> _______________________________________________
>>>> Nek5000-users mailing list
>>>> Nek5000-users at lists.mcs.anl.gov
>>>> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
>>> _______________________________________________
>>> Nek5000-users mailing list
>>> Nek5000-users at lists.mcs.anl.gov
>>> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
>> _______________________________________________
>> Nek5000-users mailing list
>> Nek5000-users at lists.mcs.anl.gov
>> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
>> _______________________________________________
>> Nek5000-users mailing list
>> Nek5000-users at lists.mcs.anl.gov
>> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
>
> _______________________________________________
> Nek5000-users mailing list
> Nek5000-users at lists.mcs.anl.gov
> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users