[Fabf]

Dear Foamers,

a simple question for someone who is more skilled in programming:

Based on the modified adjointShapeOptimization version proposed by Chalmers [1] I want to change the definition of the user defined target-outflow-direction from cartesian coordinates to a normal-to-boundary definition. Aim is to use it for a radial-diffusor configuration.

In fvSolution that user definined direction "Ud" is called from Uduserdefnodim:

objectiveFunctionDict
{
    objectiveFunction           2;
    numberObjectivePatches      1;
    objectivePatchesNames       (outlet1);
    Uduserdefnodim              (1 0 0);
}

In adjointOutletVelocityUniFvPatchVectorField.C it is used as

    const fvPatchField<vector>& Udp =
        patch().lookupPatchField<volVectorField,vector>("Ud");

    vectorField Udp_n = (Udp & patch().nf())*patch().nf();

Since I do not want the normal components from Ud as Udp_n, I just changed Udp_n to

    vectorField Udp_n = 1.0*patch().nf();

In adjointOutletPressureUni.C Ud is used as

    scalarField Udp_n = (Udp & patch().nf());

I replaced it by "1".

But that solution was obviously too easy. My case crashed after 17 iterations:

Objective function (Flow uni at outlet) J: 0
smoothSolver:  Solving for Ux, Initial residual = 0.820834, Final residual = 0.0720524, No Iterations 2
smoothSolver:  Solving for Uy, Initial residual = 0.452784, Final residual = 0.0148198, No Iterations 2
smoothSolver:  Solving for Uz, Initial residual = 0.701969, Final residual = 0.028579, No Iterations 2
GAMG:  Solving for p, Initial residual = 0.860458, Final residual = 0.00859021, No Iterations 2
time step continuity errors : sum local = 1.51941e+37, global = 9.63723e+34, cumulative = 9.63723e+34
smoothSolver:  Solving for Uax, Initial residual = 0.00135215, Final residual = 4.36927e-06, No Iterations 10
smoothSolver:  Solving for Uay, Initial residual = 0.00157789, Final residual = 3.25466e-05, No Iterations 10
smoothSolver:  Solving for Uaz, Initial residual = 0.00378074, Final residual = 1.52477e-05, No Iterations 10
GAMG:  Solving for pa, Initial residual = 0.925591, Final residual = 0.0923264, No Iterations 1000
Adjoint continuity errors : sum local = 5.07986e+17, global = -1.06656e+15, cumulative = -1.06656e+15
smoothSolver:  Solving for epsilon, Initial residual = 1, Final residual = 0.0334108, No Iterations 48
smoothSolver:  Solving for k, Initial residual = 3.60959e-14, Final residual = 3.60959e-14, No Iterations 0
ExecutionTime = 229.96 s


Time = 17

Objective function (Flow uni at outlet) J: 0
[2] #0  Foam::error::printStack(Foam::Ostream&)[11] #0  Foam::error::printStack(Foam::Ostream&)[5] #0  Foam::error::printStack(Foam::Ostream&)[1] #0  Foam::error::printStack(Foam::Ostream&)[0] #0  Foam::error::printStack(Foam::Ostream&)[9] #0  Foam::error::printStack(Foam::Ostream&) at ??:?
[5] #1  Foam::sigFpe::sigHandler(int) at ??:?
[1] #1  Foam::sigFpe::sigHandler(int) at ??:?
[11] #1  Foam::sigFpe::sigHandler(int) at ??:?
[0] #1  Foam::sigFpe::sigHandler(int) at ??:?
[2] #1  Foam::sigFpe::sigHandler(int) at ??:?
[9] #1  Foam::sigFpe::sigHandler(int) at ??:?
[1] #2  
 at ??:?
[11] #2  
 at ??:?
[0] #2  
 at ??:?
[5] #2  
 at ??:?
[2] #2  
 at ??:?
[9] #2  
[1]  at sigaction.c:?
[1] #3  Foam::GaussSeidelSmoother::smooth(Foam::word const&, Foam::Field<double>&, Foam::lduMatrix const&, Foam::Field<double> const&, Foam::FieldField<Foam::Field, double> const&, Foam::UPtrList<Foam::lduInterfaceField const> const&, unsigned char, int)[11]  at sigaction.c:?
[11] #3  [0]  at sigaction.c:?
[0] #3  Foam::GaussSeidelSmoother::smooth(Foam::word const&, Foam::Field<double>&, Foam::lduMatrix const&, Foam::Field<double> const&, Foam::FieldField<Foam::Field, double> const&, Foam::UPtrList<Foam::lduInterfaceField const> const&, unsigned char, int)Foam::GaussSeidelSmoother::smooth(Foam::word const&, Foam::Field<double>&, Foam::lduMatrix const&, Foam::Field<double> const&, Foam::FieldField<Foam::Field, double> const&, Foam::UPtrList<Foam::lduInterfaceField const> const&, unsigned char, int)[2]  at sigaction.c:?
[2] #3  Foam::GaussSeidelSmoother::smooth(Foam::word const&, Foam::Field<double>&, Foam::lduMatrix const&, Foam::Field<double> const&, Foam::FieldField<Foam::Field, double> const&, Foam::UPtrList<Foam::lduInterfaceField const> const&, unsigned char, int)[5]  at sigaction.c:?
[5] #3  Foam::GaussSeidelSmoother::smooth(Foam::word const&, Foam::Field<double>&, Foam::lduMatrix const&, Foam::Field<double> const&, Foam::FieldField<Foam::Field, double> const&, Foam::UPtrList<Foam::lduInterfaceField const> const&, unsigned char, int)[9]  at sigaction.c:?
[9] #3  Foam::GaussSeidelSmoother::smooth(Foam::word const&, Foam::Field<double>&, Foam::lduMatrix const&, Foam::Field<double> const&, Foam::FieldField<Foam::Field, double> const&, Foam::UPtrList<Foam::lduInterfaceField const> const&, unsigned char, int) at ??:?
[1] #4  Foam::GaussSeidelSmoother::smooth(Foam::Field<double>&, Foam::Field<double> const&, unsigned char, int) const at ??:?
[0] #4  Foam::GaussSeidelSmoother::smooth(Foam::Field<double>&, Foam::Field<double> const&, unsigned char, int) const at ??:?
[11] #4  Foam::GaussSeidelSmoother::smooth(Foam::Field<double>&, Foam::Field<double> const&, unsigned char, int) const at ??:?
[2] #4  Foam::GaussSeidelSmoother::smooth(Foam::Field<double>&, Foam::Field<double> const&, unsigned char, int) const at ??:?
[9] #4  Foam::GaussSeidelSmoother::smooth(Foam::Field<double>&, Foam::Field<double> const&, unsigned char, int) const at ??:?
[5] #4  Foam::GaussSeidelSmoother::smooth(Foam::Field<double>&, Foam::Field<double> const&, unsigned char, int) const at ??:?
[1] #5  Foam::smoothSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
[11] #5  Foam::smoothSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
[0] #5  Foam::smoothSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
[9] #5  Foam::smoothSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
[2] #5  Foam::smoothSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
[5] #5  Foam::smoothSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
[1] #6   at ??:?
[11] #6  
 at ??:?
[0] #6   at ??:?
[9] #6   at ??:?
[5] #6   at ??:?
[2] #6  [1]  at ??:?
[1] #7  




[11]  at ??:?
[11] #7  
[0]  at ??:?
[0] #7  [5]  at ??:?
[5] #7  [1]  at ??:?
[1] #8  [9]  at ??:?
[9] #7  [2]  at ??:?
[2] #7  




[0]  at ??:?
[0] #8  [5]  at ??:?
[5] #8  
[1]  at ??:?
[1] #9  [2]  at ??:?
[2] #8  [11]  at ??:?
[11] #8  

[9]  at ??:?
[9] #8  


[0]  at ??:?
[0] #9  [5]  at ??:?
[5] #9  
[2]  at ??:?
[2] #9  [1]  at ??:?
[1] #10  __libc_start_main[11]  at ??:?
[11] #9  

[9]  at ??:?
[9] #9   at ??:?
[1] #11  
[5]  at ??:?
[5] #10  __libc_start_main
[2]  at ??:?
[2] #10  __libc_start_main
[0]  at ??:?
[0] #10  __libc_start_main
[11]  at ??:?
[11] #10  __libc_start_main at ??:?
[5] #11   at ??:?
[2] #11  [1]  at ??:?
[thermo-04:22271] *** Process received signal ***
[thermo-04:22271] Signal: Floating point exception (8)
[thermo-04:22271] Signal code:  (-6)
[thermo-04:22271] Failing at address: 0xdc92000056ff
[9]  at ??:?
[9] #10  __libc_start_main at ??:?
[0] #11  [thermo-04:22271] [ 0] /lib64/libc.so.6() [0x3c3e2326a0]
[thermo-04:22271] [ 1] /lib64/libc.so.6(gsignal+0x35) [0x3c3e232625]
[thermo-04:22271] [ 2] /lib64/libc.so.6() [0x3c3e2326a0]
[thermo-04:22271] [ 3] /opt/OpenFOAM/OpenFOAM-2.3.0/platforms/linux64Gcc48DPOpt/lib/libOpenFOAM.so(_ZN4Foam19GaussSeidelSmoother6smoothERKNS_4wordERNS_5FieldIdEERKNS_9lduMatrixERKS5_RKNS_10FieldFieldIS4_dEERKNS_8UPtrListIKNS_17lduInterfaceFieldEEEhi+0x345) [0x7f4483ddbca5]
[thermo-04:22271] [ 4] /opt/OpenFOAM/OpenFOAM-2.3.0/platforms/linux64Gcc48DPOpt/lib/libOpenFOAM.so(_ZNK4Foam19GaussSeidelSmoother6smoothERNS_5FieldIdEERKS2_hi+0x2d) [0x7f4483ddbf0d]
[thermo-04:22271] [ 5] /opt/OpenFOAM/OpenFOAM-2.3.0/platforms/linux64Gcc48DPOpt/lib/libOpenFOAM.so(_ZNK4Foam12smoothSolver5solveERNS_5FieldIdEERKS2_h+0x4f7) [0x7f4483dd5937]
[thermo-04:22271] [ 6] my_UNAdjointShapeOptimizationFoam() [0x459a4a]
[thermo-04:22271] [ 7] my_UNAdjointShapeOptimizationFoam() [0x46ce21]
[thermo-04:22271] [ 8] my_UNAdjointShapeOptimizationFoam() [0x46d0f6]
[thermo-04:22271] [ 9] my_UNAdjointShapeOptimizationFoam() [0x428d30]
[thermo-04:22271] [10] /lib64/libc.so.6(__libc_start_main+0xfd) [0x3c3e21ed5d]
[thermo-04:22271] [11] my_UNAdjointShapeOptimizationFoam() [0x42bb41]
[thermo-04:22271] *** End of error message ***

Has somebody a better solution for this? Do I have to change more than this tiny bit of code? I know, that there is bc called surfaceNormalFixedValue, which basically makes what I want. I looked at the code but I don't know how to use it exactly for my problem...

Thanks in advance!


[1] http://www.tfd.chalmers.se/~hani/kur...ortAdjoint.pdf

[Bloerb]

The version of this solver supplied with an standard openfoam installation works just fine. From my point of view the chalmers version is wrongly implemented in some minor and major points.

Take a look at equation 2.2 in your linked paper and the missing part in the mentioned papers. The standard solver optimizes for total pressure loss (p+1/2v²) between inlet and outlet. The missing terms in adjointOutletVelocity and adjointOutletPressure are left out because they are of second order. Not only is the viscosity a usually small value, but those gradients of the adjoint near the outlet are negligible. Hence the optimization goal is the same, but these small values are neglected for a more robust boundary condition. Adding these terms as done in the chalmers paper is possible, but you should do some calculations and comparisons if it is necessary. It lowers your stability and forces you to use robin type boundaries. Nevertheless chalmers implemented those changes correctly.

The chalmers paper changes the sign of the gradient step. This is wrong and only works in this paper because the inflow boundary condition is implemented wrongly. Note that in the tutorial case of openfoam the inlet boundary is (-1 0 0). The inlet of the normal system is the outlet of the adjoint. The adjoint system flows backwards with supplied velocity at the "outlet". This is not done in the chalmers papers and wrong. They use u_n=v_n while it should be u_n=- v_n

The geometry is changed based on u*v. A scalar product of adjoint and normal velocity. In 2d this means the angle between them is the only thing that determines the sign of the sensitivity. negative for removing or positive for adding material. Please take an indepth look at the difference between the way chalmers does this and how the mathematics of othmer describes it. You can see that changing the sign of the gradient step counteracts the boundary. Nevertheless this is not correct imo.

My recommendation to you is changing the direction of the inflow boundary. Afterwards change the step size lambda to zero and look at the solution of the adjoint. I'd even initialize it with a converged solution of the primary system. You should be able to tell why your solver blows uo. This way you are only solving the adjoint system and should look at how it flows and how its pressure and velocity values change over the iterations. On the first look your change looks decent to me. I am tired though

[Fabf]

Hello Bloerb,

thank you very much for your detailed and very helpfull explanation!

To 1) All right, this probably explaines why I recently had issues with full-blocking of the field or other stability problems with the chalmers solver for optimizing pressure loss, while the standard implementation from othmer worked just fine. So maybe I just delete those 2nd order terms in Chalmers version. I am gonna have a closer look at the mathematics.

To 2) Oh yes. My case was based on the stanard tutorial case, so my adjoint inflow direction was outwards- I changed the sign and the results seem to be reasonable! So I am gonna just change the sign in my chalmers-based-version as well. In the chalmers-documentation they wrote about some sign-mistakes but also said, it wouldn't have any big influence --> wrong!

To 3) Okay, my results so far are not that bad. But I really have to go deeper in both implementations to understand the differences and if my results are trustable.

4) for my first mesh (blockMesh 2x2mm, 2D axisymmetric) I did not have to set lambda to zero and also a converged solution for initializing was not necessary. I just had to set alpha=1e6. Nevertheless, with my fine mesh (0,5mm) I do get total different results with those settings. Now I will try with your suggestions.

Seems to me that you are quiet experienced with this topic. Do you have some publications, github or whatever where you share your knowledge resp. your own solvermodifications?

[Bloerb]

To start you should look at othmers papers or the following by hinterberger and oelsen. "Industrial application of continuous adjoint Flow solver for the optimization of automotive exhaust systems." They give you a derivation of the adjoint equations used. Combine this with a good book about optimization, take a look at the method of lagrange multipliers and you should understand the math behind the method.

Now to the solver. Your primary solver is identical to simpleFoam. The only difference is the added darcy term used to block your velocities. Your adjoint solver is hence also similar but has the characteristics of an adjoint solver. (for example reversed convective transport) adjointshapeoptimizationfoam now does the follwoing

calll primary solver to solve for u and p
call adjoint solver for ua und pa
use u and ua to update the design parameter field alpha.

Those are essentially 3 solvers. Lets say primaryFoam, adjointFoam and topologyUpdateFoam. In this implementation each is solved for one iteration. You could also fully converge both calculations before changing the topology. In fact this would be the mathematical way. It helps to understand the adjoint solution. It is however often not that easy to get a converging adjoint solution.

This is important because there are three steps for good results. Check if your primary solution is decent. If it is check if your adjoint solution behaves as you'd expect. Now changing the topology in a meaningful way is only possible if both results are plausible. Hence i'd set lambda to zero until you are certain both solvers do what they are supposed to. Add the sensitivity output like chalmers did and check if those look plausible. And do this for the easiest geometry you can think of. One where your goal is easy to check, like the tutorial case provided. Once your sensitivity is what it should be, change the topology. This is a topic on its own.

As for alpha=1e6. Because the update of the topology is done by u*ua you might see bad results with such a high value. For an infinite value of alpha there is no velocity. u=0. and also ua=0. Hence you can not lower your value of alpha. Your sensitivity is zero. In a nutshell: It is much easier to add solid regions. Removing them becomes nearly impossible. A low value on the other hand might give you soft walls. The difference between solid and fluid becomes a gray mess instead of a black and white solution. Logically because you still have velocity left in your solid regions. Setting this, the relaxation factor, the step size, etc are all trial and error.

Sadly there is no code i can share. I have added several things to this solver, but no finished version yet. But if you have any questions don't hesitate to contact me. I am still learning about this myself though and hence no expert. But we might be able to help each other.

[Fabf]

I had some time to dig through the literature and do some further steps regarding my case.

I took the original code from othmer and added
- the sensitivity output
- the cost function printing during run-time (was implemented wrong by chalmers, but described correctly in the documentation). Anyways, works not in parallel.
- added the uniformity-BCs in a similiar way like chalmers did. Created the solver versions: with and without 2nd-order terms (called "robust Solver")
- changed the uniformity-BCs according to normal outflow w.rsp.t. outlet boundary patch.

Starting from a simple testcase (pitzdaily) everthing worked fine. There was no difference between the robust solver and the normal one.

Happy with that result I switched to my case (something like a radial diffusor, simulated as 2D-axisymmetric). After initializing the primal solution I run the adjoint solution and updated alpha.

Now the results are quite different to the testcase. In most of the domain the solution looks really good for the first few thousend steps (lambda=1e1..1e5). Then the domain with low speed gets blocked, just a few channels remain. What confuses me: Near the outlet patch the blockage is growing like a ramp, like it wants to have a higher velocity there. But the average velocity is much higher than what I defined. Also the flow is not uniform anymore. I can see an adjoint velocity in those areas, so its an issue with the BC-implementation? I know, that the adjoint solution is "growing" from the patches into the domain. Do you have any explanation for that behaviour? Besides: As for the testcase, there was no big difference between the robust and the normal solver (just in detail some cells were blocked, where others not).

Another issue: Since the k-epsilon-modell gives wrong results for that geometry (stagnation flow) I used k-omega SST model for the primal solution. Worked fine so far. At this point I recognized, that I did not really understood the concept of frozen turbulence and what influence it has on my solution. Since the turbulence parameters are computed and also used at some points by the adjoint solver (see 2nd Order terms) - what is the sense behind "freezing"?

I played around a bit with the schemes. Although Upwind is recommended I had no problems by using limitedLinear or LUST. So it's worth a try.

See the attached file for the porosity field. Outlet is right. Inlet is at the bottom. Axis of symmetry is left.

[verboomj]

Hi Fabf and Bloerb,

Thanks for this in depth discussion, It has been very useful for my understanding of the implementation of this method.

In addition to the papers mentioned here I've found and extension to the current frozen turbulence assumption, a 2016 paper. The paper is quite thorough on the objectives and the adjoint equations. In chapter 7 of the paper they apply the method on a nice 2D problem which is rather simple to use as a benchmark. paper: http://link.springer.com/article/10....831-014-9141-9

As far as I understand the turbulence is frozen to simplify the model, but this is only done for the adjoint equations. The frozen turbulence assumption essentially says that the influence to the convection by the turbulence can be neglected because it is small w.r.t. the mean advection speed. It holds if the turbulence intensity is rather small.. http://glossary.ametsoc.org/wiki/Taylor's_hypothesis

Since this assumption is used for the adjoint equations, your sensitivities found by the adjoint method are not exact, however it is a good approximation as long as the turbulence intensity is not to large.

I've currently started my master thesis project and will be working with this solver in OpenFOAM and finally apply it to a rather simple 3D problem.
Anyway, I've done the derivations by Othmer and looked at the implementations so far. As I will be busy on this topic for the coming 8 months it might be nice to have some fellow FOAMERS busy with this solver.

If some questions arise, and you guys don't mind it might be nice to discuss some topics with each other.

[palarcon]

Is anybody working on that know?
I'm a PhD student, specialized on topology optimization, but OpenFOAM is not my strong point and I want to perform topology optimization using OpenFOAM, mainly improve the optimization part of the OpenFOAM solver.
If anybody is interested please leave a message.