[openfoam_aero]

Dear all,

OpenFOAM version: OpenFOAM 7

I am currently working on data assimilation and I have to solve a set of adjoint equations for a domain consisting of a flow past a cylinder. The thing is, these equations are similar to NS equations but with the exception that they are not non-linear. I present the equattions below
adjoint_equation.png

In these equations, I have the value of \Tilde{u} and the term on the RHS (which is a source term in this case). I need to solve for u_dagger and p_dagger. There is a continuity equation that is also present given by
continuity.png

Now I use the icoFoam solver by removing the time derivative (since it is steady state problem). I am presenting the code below

#include "fvCFD.H"
#include "simpleControl.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    #include "setRootCaseLists.H"
    #include "createTime.H"
    #include "createMesh.H"

    simpleControl simple(mesh);

    #include "createFields.H"
    #include "initContinuityErrs.H"

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    Info<< "\nStarting time loop\n" << endl;

    while (simple.loop(runTime))
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        //#include "CourantNo.H"

        // Momentum predictor


	tmp<fvVectorMatrix> tUEqn
	(
		 - (fvm::laplacian(nu,UDagger)) +(U1 & T(fvc::grad(UDagger))) - (UDagger & fvc::grad(U1)) - dEdU
	
	);

        
        fvVectorMatrix& UEqn = tUEqn.ref();

    	UEqn.relax();


        solve(UEqn == -fvc::grad(pDagger));
        


        volScalarField rAU(1.0/UEqn.A());
        
        volVectorField HbyA(constrainHbyA(rAU*UEqn.H(),UDagger, pDagger));
        
        surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
	
	adjustPhi(phiHbyA, UDagger, pDagger);

	
	tmp<volScalarField> rAtU(rAU);

    	if (simple.consistent())
    		{
		rAtU = 1.0/(1.0/rAU - UEqn.H1());
		phiHbyA +=
		    fvc::interpolate(rAtU() - rAU)*fvc::snGrad(pDagger)*mesh.magSf();
		HbyA -= (rAU - rAtU())*fvc::grad(pDagger);
    		}
        
        tUEqn.clear();
        


        // Update the pressure BCs to ensure flux consistency
        constrainPressure(pDagger, UDagger, phiHbyA, rAU);

        // Non-orthogonal pressure corrector loop
        //while (simple.correctNonOrthogonal())
        //{
             // Pressure corrector

             fvScalarMatrix pEqn
             (
                 fvm::laplacian(rAU, pDagger) == fvc::div(phiHbyA)
             );

 //            pEqn.setReference(pRefCell, pRefValue);

             pEqn.solve();

             //if (simple.finalNonOrthogonalIter())
             //{
                 phi = phiHbyA - pEqn.flux();
             //}
         //}

         #include "continuityErrs.H"
            
         pDagger.relax();

         UDagger = HbyA - rAU*fvc::grad(pDagger);
         UDagger.correctBoundaryConditions();

         runTime.write();

         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }

    Info<< "End\n" << endl;

    return 0;
}

The boundary conditions for U_Dagger and p_dagger are similar to the flow past a cylinder case with symmetry BCs on the top and bottom, empty BCs for front and back, noSlip for the cylinder, inlet BC for velocity and zeroGradient outlet for velocity.

I use a GAMG solver for p_dagger and PBiCG solver for U_Dagger with residuals of 10^-4 for each.

I get a U_Dagger field after about 600 iterations but it is not correct when compared to the reference case from literature. The gradients close to the cylinder wall are not present in my case. I tried debugging a lot and looked through a lot of options.

My questions
1. Am I setting up the equations correctly?
2. Am I using the right method of solving these equations? Is there a simpler alternative?

I suspect that the issue might be with fvVectorMatrix tUEqn. Any suggestions on the problem will be appreciated.
P.S i have included #createPhi1.H at the end of my #createFields.H file and the #createPhi1.H file looks like this

Info<< "Reading/calculating face flux field phi\n" << endl;

surfaceScalarField phi
(
    IOobject
    (
        "phi",
        runTime.timeName(),
        mesh,
        IOobject::READ_IF_PRESENT,
        IOobject::AUTO_WRITE
    ),
    fvc::flux(UDagger)
);