CFD Online Discussion Forums - Volume integration of solution variable

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- SU2 (https://www.cfd-online.com/Forums/su2/)

- - Volume integration of solution variable (https://www.cfd-online.com/Forums/su2/166391-volume-integration-solution-variable.html)

Volume integration of solution variable

Hi!

I need to compute a volume integral over the computational domain (for now a 2D airfoil). Specifically, I need to find {\int \int - x \nabla\cdot V dx dy}.

Now I could not find any built-in capability to do so, so I have to implement it myself.

Could you please give some advice on best way of doing that? My guess is that I can reuse some parts of the current code where the residuals are being computed.

Thanks!

Hi!

In this case my recommendation is to implement the formula at the level of the solver class. For example, dealing with Euler equations, you can use

void CEulerSolver::Postprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh) { } in solver_direct_mean.cpp that it is empty.

To find the name of the variables, it is useful to take a look at other well known subroutines like

void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_container, CNumerics *numerics, CConfig *config, unsigned short iMesh).

In your particular case you will need something like

void CEulerSolver::Postprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config,
unsigned short iMesh) {

unsigned long iPoint;
su2double Area, dvx_dx, dvx_dy, dvy_dx, dvy_dy, x;

for (iPoint = 0; iPoint < nPoint; iPoint ++) {
Area = geometry->node[iPoint]->GetVolume();
x = geometry->node[iPoint]->GetCoord(0);
dvx_dx = node[iPoint]->GetGradient_Primitive()[1][0];
dvx_dy = node[iPoint]->GetGradient_Primitive()[1][1];
dvy_dx = node[iPoint]->GetGradient_Primitive()[2][0];
dvy_dy = node[iPoint]->GetGradient_Primitive()[2][1];

// The formula should be computed here
// a simple cout will output the result

}

}

Thanks for using SU2,
Best,
Francisco

Hi Francisco,

Thanks for your answer. However, if I use this, it computes it every MG iteration, so I moved it to output_structure.cpp.

Regardless of that, the gradients always return zero. Where are the velocity gradients actually computed/which function should I call for them?

I wrote the function

Code:

 void COutput::CompDoubletCirc(CGeometry ***geometry,

                CSolver ****solver_container,

                CConfig **config,

                unsigned short val_iZone, unsigned long iExtIter)

where the main loop is this (note that kappax is not computed correctly now, but I want to see whether any of the derivatives is nonzero -- I have also tested them separately):

Code:

for (iPoint = 0; iPoint < geometry[val_iZone][FinestMesh]->GetnPoint(); iPoint ++) {



                        Area         = geometry[val_iZone][FinestMesh]->node[iPoint]->GetVolume();

                        x                 = geometry[val_iZone][FinestMesh]->node[iPoint]->GetCoord(0);



                        dvx_dx = solver_container[val_iZone][FinestMesh][FLOW_SOL]->node[iPoint]->GetGradient_Primitive()[1][0];

                        dvx_dy = solver_container[val_iZone][FinestMesh][FLOW_SOL]->node[iPoint]->GetGradient_Primitive()[1][1];

                        dvy_dx = solver_container[val_iZone][FinestMesh][FLOW_SOL]->node[iPoint]->GetGradient_Primitive()[2][0];

                        dvy_dy = solver_container[val_iZone][FinestMesh][FLOW_SOL]->node[iPoint]->GetGradient_Primitive()[2][1];



                        // Compute doublet strength (method 1)

                        //                 \kappa_x = \int \int x \nabla \cdot V dA



                        kappax = kappax + dvx_dx + dvx_dy + dvy_dx + dvy_dy; // Area*x*(dvx_dx + dvy_dy);



                }

And I call this function in SU2_CFD.cpp with

Code:

 output->CompDoubletCirc(geometry_container, solver_container,

                  config_container, ZONE_0, ExtIter);

Thanks!

I see...

Please notice that the gradients are computed in

void CEulerSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh, unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output) {

using

SetPrimitive_Gradient_GG(geometry, config);

or

SetPrimitive_Gradient_LS(geometry, config);

If you look at the code, SU2 doesn't compute gradients for centered schemes like JST. Are you running JST?

If you want to be sure that you have computed the gradients maybe you can call the subroutine

solver_container[val_iZone][FinestMesh][FLOW_SOL]->SetPrimitive_Gradient_GG(geometry, config);

just before running your code.

Best,
Francisco

Awesome, that was exactly the problem, I was using JST and the gradients weren't being computed.

Thanks!