CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

 Envolly June 19, 2013 16:39

2 Attachment(s)
Hello I am new here and this is my first post. I am a mechanical engineering student and currently writing my master’s thesis. This is about the adjoint optimization method.

I want to change the optimization target form the current implementation in OF 2.2.0 (minimize pressure-losses) to the target of the flow uniformity at the outlet.
Therefore I have tried to implement the equations for this cost function based on the papers of Othmer and Hinterberger/Olesen.

A continuous adjoint formulation for the computation of topological and surface sensitivities of ducted flows (C. Othmer)

Automatic Geometry Optimization of Exhaust Systems Based on Sensitivities Computed by a Continuous Adjoint CFD Method in OpenFOAM (Hinterberger, C. and Olesen, M.)

According to Othmer only the code of adjointOutletPressure and adjointOutletVelocity has to be the changed. I have done this as follows:

Code:

``` const fvsPatchField<scalar>& phiap =         patch().lookupPatchField<surfaceScalarField, scalar>("phia");     const fvPatchField<vector>& Up =         patch().lookupPatchField<volVectorField, vector>("U");     const fvPatchField<vector>& Udp =         patch().lookupPatchField<volVectorField, vector>("Ud");     scalarField Un(mag(patch().nf() & Up));            % normal component primal vel.     scalarField Udn(mag(patch().nf() & Udp));          % normal component desired vel     vectorField Ut(Up - patch().nf()*Un);              % tangential component primal vel.     vectorField Udt(Udp - patch().nf()*Udn);            % tangential component desired vel.     vectorField UtHat((Ut - Udt)/(Un + SMALL));        % tangential component adjoint vel.     vectorField Ua(phiap*patch().Sf()/sqr(patch().magSf()));  % normal component adjoint vel.     vectorField::operator=(Ua - UtHat);```

Code:

```const fvsPatchField<scalar>& phip =         patch().lookupPatchField<surfaceScalarField, scalar>("phi");     const fvsPatchField<scalar>& phiap =         patch().lookupPatchField<surfaceScalarField, scalar>("phia");     const fvPatchField<vector>& Up =         patch().lookupPatchField<volVectorField, vector>("U");     const fvPatchField<vector>& Uap =         patch().lookupPatchField<volVectorField, vector>("Ua");     const fvPatchField<vector>& Udp =         patch().lookupPatchField<volVectorField, vector>("Ud");     scalarField Udn(mag(patch().nf() & Udp));     operator==((phiap/patch().magSf() + 1.0)*phip/patch().magSf() + (Up & Uap) - Udn);```
The vector Field Ud contains the boundary values of the desired velocity. U is the primal velocity and Ua the adjoint velocity.

I have applied this code to a simple (turbulent) 2d example with one inlet and three outlets. The primal velocity at the inlet was set to 1 m/s and the desired velocity at the outlet to 0.667 m/s. The boundaries for the primal variables and the adjoint-pressure were set as in the tutorial case. The adjoint-velocity at the inlet was set to zero as stated in the paper of Othmer.
The convergence was quite good but the result is strange as you can see in the picture "3outlets".

I have tried another simple example of a laminar 2d duct flow. The primal velocity and the adjoint velocity were chosen to be equal. Normally, no cells should be blocked in that case.
The result is that nearly all cells are blocked. Just a small flow passage remains (picture: duct)

My questions now are:

- is the implementation of the equations for the optimization target right?
- do I have to change something else in the code ?
- according to Othmer, the difference of the primal and desired velocity would act as a source term in the outlet bc. How can this be understood?
- any other notes ?

best regards

 sylvester June 20, 2013 08:55

Quote:
 Originally Posted by Envolly (Post 434901) - is the implementation of the equations for the optimization target right? - do I have to change something else in the code ? - according to Othmer, the difference of the primal and desired velocity would act as a source term in the outlet bc. How can this be understood?
Hi Envolly,

Welcome to the forum and the world of adjoint optimization!

I have experimented with this a bit before, but I never have found the time to fully test my code, so you should probably treat the following with care.

1. My bc's are comparable to yours, although some parts I did differently. For example, I simplified the velocity boundary condition by assuming only normal desired and adjoint flow. Also, looking at the pressure bc I have +Udn, where you have -Udn. (But looking at my code again, I am pretty sure the error is mine!)

2. Basically, only the bc's need to be adjusted to fit your cost function. However, for stability often other modifications are necessary; see some other posts by me about this subject. Stability is, still, an ongoing problem for me, therefore I cannot help you more on this than with what I have written before.

3. The outlet bc acts as the adjoint source as the propagation of the adjoint variables is reversed with respect to your primal flow.

I hope this helps you.

Good luck!

regards,
Sylvester

 Envolly July 1, 2013 03:01

Hi Sylvester,

To 1:
Quote:
 My bc's are comparable to yours, although some parts I did differently. For example, I simplified the velocity boundary condition by assuming only normal desired and adjoint flow. Also, looking at the pressure bc I have +Udn, where you have -Udn. (But looking at my code again, I am pretty sure the error is mine!)
According to the papers I have quoted above, it should be minus the normal component of the desired velocity (-Udn).
I have chosen a desired velocity which has only a normal component, too. If I simply
the bc for the adjoint velocity it look like:

Code:

` vectorField ::operator=(Ua)`
The pressure boundary stays unchanged.
Unfortunately I get similar strange results like before.

To 2:
Quote:
 Basically, only the bc's need to be adjusted to fit your cost function.
This agrees with the statements in the articles I have quoted above and therefore my code should be ok.

Quote:
 However, for stability often other modifications are necessary; see some other posts by me about this subject. Stability is, still, an ongoing problem for me, therefore I cannot help you more on this than with what I have written before.

1. To start from a converged primal solution
2. Update the alpha-field after the adjoint field has roughly converged
3. Mesh quality
4. Pseudo staggered approach

To

3: The mesh quality in my case is not the problem
4: Is that much programming effort ?

1+2: I have already tried to start from a converged primal solution without a positive effect.

I will try to start the adjoint method on a converged primal solution without updating the alpha field to get a converged adjoint field.

However, when I look at the boundaries I have some doubts how the adjoint flow field should arise or look like (see attachment).
The conditions for Ua at the inlet is (accoriding to Othmer) a fixedValue (Ua =0 ) and the the one at the outlet
acts as a fixedValue, or ?

The vector plot (3outlets) of my first case shows what I mean. The adjoint flow enters the domain at the first outlet and leaves it trough the others. Could, or should that be right ?

Best regards

Envolly

 immortality July 2, 2013 03:55

hi
I'm unfamiliar with adjoint method but a question occurred to me.
Are you implementing a fixedValue of velocity on the output instead of fixedValue pressure?
Does it seem right?

 Envolly July 2, 2013 05:49

Hi immortality,

as far as I could determine the type of boundaries for the adjoint variables ( my c++ knowledge is not so good), both are fixedValues.

You will find for the pressure:

Code:

```    operator==((phiap/patch().magSf() - 1.0)*phip/patch().magSf() + (Up & Uap));     fixedValueFvPatchScalarField::updateCoeffs();```
and for the velocity:
Code:

```vectorField::operator=(phiap*patch().Sf()/sqr(patch().magSf()) + UtHat);     //vectorField::operator=(Uan + UtHat);     fixedValueFvPatchVectorField::updateCoeffs();```
The only difference from my point is the way how the values are assigned.

pressure:

Code:

`operator==()`
velocity:

Code:

`vectorField::operator=()`
Is there a difference in processing the values​​.

Best regards

Envolly

 immortality July 2, 2013 06:21

then U and p are set on outlet.why it has be done in this manner?
and what means that formula in mathematical form?where they come from?do you know?

 Envolly July 3, 2013 04:28

For detailed description I would refer you to the papers which I have quoted in my first post:

Quote:
 A continuous adjoint formulation for the computation of topological and surface sensitivities of ducted flows (C. Othmer) Automatic Geometry Optimization of Exhaust Systems Based on Sensitivities Computed by a Continuous Adjoint CFD Method in OpenFOAM (Hinterberger, C. and Olesen, M.)

Ua and Pa (adjoint vel. and adjoint pressure) are the Lagrange multipliers.
The adjoint-method uses a Lagrange function because it is a constrained optimization problem.

The optimization target (total pressure losses, energy dissipation, flow uniformity)
is defined in the boundary conditions. That´s why they are so important.

regards

Envolly

 immortality July 3, 2013 08:34

so they aren't real pressure and velocity,right?

 Envolly July 9, 2013 01:29

from my point of view, they have no real physical meaning. the adjoint velocity is needed to calculate the sensitivity.

But if I'm wrong please anybody correct me..

 All times are GMT -4. The time now is 09:33.