[Envolly]

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:

adjointOutletVelocity:

 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);

adjointOutletPressure:

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 ?


Thanks a lot in advance, appreciate all answers

best regards

[sylvester]

Originally Posted by Envolly

- 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]

Hi Sylvester,

thanks für your reply.

To 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!)

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:

 vectorField ::operator=(Ua)

The pressure boundary stays unchanged.
Unfortunately I get similar strange results like before.

To 2:

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.

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.

I have read some of your posts before. Some of your instructions are:

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]

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]

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.
Have a look in the adjointOutletPressure.C and adjointOutletVelocity.C.

You will find for the pressure:

    operator==((phiap/patch().magSf() - 1.0)*phip/patch().magSf() + (Up & Uap));

    fixedValueFvPatchScalarField::updateCoeffs();

and for the velocity:

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:

operator==()

velocity:

vectorField::operator=()

Is there a difference in processing the values​​.


Best regards

Envolly

[immortality]

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]

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

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.)

short answer:

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]

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

[Envolly]

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..