|
[Sponsors] |
February 8, 2016, 16:01 |
Volume integration of solution variable
|
#1 |
New Member
Anonymous
Join Date: Aug 2013
Location: Europe
Posts: 24
Rep Power: 13 |
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! |
|
February 13, 2016, 11:05 |
|
#2 |
Super Moderator
Francisco Palacios
Join Date: Jan 2013
Location: Long Beach, CA
Posts: 404
Rep Power: 15 |
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 |
|
February 28, 2016, 20:39 |
|
#3 |
New Member
Anonymous
Join Date: Aug 2013
Location: Europe
Posts: 24
Rep Power: 13 |
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) 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); } Code:
output->CompDoubletCirc(geometry_container, solver_container, config_container, ZONE_0, ExtIter); |
|
February 28, 2016, 21:08 |
|
#4 |
Super Moderator
Francisco Palacios
Join Date: Jan 2013
Location: Long Beach, CA
Posts: 404
Rep Power: 15 |
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 |
|
February 29, 2016, 00:40 |
|
#5 |
New Member
Anonymous
Join Date: Aug 2013
Location: Europe
Posts: 24
Rep Power: 13 |
Awesome, that was exactly the problem, I was using JST and the gradients weren't being computed.
Thanks! |
|
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
alphaEqn.H in twoPhaseEulerFoam | cheng1988sjtu | OpenFOAM Bugs | 15 | May 1, 2016 17:12 |
[blockMesh] non-orthogonal faces and incorrect orientation? | nennbs | OpenFOAM Meshing & Mesh Conversion | 7 | April 17, 2013 06:42 |
On the damBreak4phaseFine cases | paean | OpenFOAM Running, Solving & CFD | 0 | November 14, 2008 22:14 |
UDF solution variable for time-averaged mass frac? | A. S. | FLUENT | 0 | May 14, 2007 17:44 |
Env variable not set | gruber2 | OpenFOAM Installation | 5 | December 30, 2005 05:27 |