# Solving a coupled linear PDE using modified icoFoam solver

 User Name Remember Me Password
 Register Blogs Members List Search Today's Posts Mark Forums Read

 June 7, 2022, 04:12 Solving a coupled linear PDE using modified icoFoam solver #1 Member   Uttam Join Date: May 2020 Location: Southampton, United Kingdom Posts: 35 Rep Power: 6 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 Code: ```#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 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 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 Code: ```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) );```

 Tags icofoam, linear pde, simple algorithm

 Thread Tools Search this Thread Search this Thread: Advanced Search Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are Off Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post jyotir OpenFOAM Running, Solving & CFD 0 November 17, 2021 10:49 Herwig OpenFOAM Running, Solving & CFD 17 November 19, 2019 13:47 lpz_michele OpenFOAM Running, Solving & CFD 53 October 19, 2015 02:50 nileshjrane OpenFOAM Running, Solving & CFD 4 February 25, 2013 04:13 liugx212 OpenFOAM Running, Solving & CFD 3 January 4, 2006 18:07

All times are GMT -4. The time now is 22:44.

 Contact Us - CFD Online - Privacy Statement - Top