CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > SU2

Volume integration of solution variable

Register Blogs Community New Posts Updated Threads Search

Like Tree1Likes
  • 1 Post By fpalacios

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   February 8, 2016, 16:01
Default Volume integration of solution variable
  #1
New Member
 
Anonymous
Join Date: Aug 2013
Location: Europe
Posts: 24
Rep Power: 13
maero21 is on a distinguished road
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!
maero21 is offline   Reply With Quote

Old   February 13, 2016, 11:05
Default
  #2
Super Moderator
 
Francisco Palacios
Join Date: Jan 2013
Location: Long Beach, CA
Posts: 404
Rep Power: 15
fpalacios is on a distinguished road
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
fpalacios is offline   Reply With Quote

Old   February 28, 2016, 20:39
Default
  #3
New Member
 
Anonymous
Join Date: Aug 2013
Location: Europe
Posts: 24
Rep Power: 13
maero21 is on a distinguished road
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!
maero21 is offline   Reply With Quote

Old   February 28, 2016, 21:08
Default
  #4
Super Moderator
 
Francisco Palacios
Join Date: Jan 2013
Location: Long Beach, CA
Posts: 404
Rep Power: 15
fpalacios is on a distinguished road
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
maero21 likes this.
fpalacios is offline   Reply With Quote

Old   February 29, 2016, 00:40
Default
  #5
New Member
 
Anonymous
Join Date: Aug 2013
Location: Europe
Posts: 24
Rep Power: 13
maero21 is on a distinguished road
Awesome, that was exactly the problem, I was using JST and the gradients weren't being computed.

Thanks!
maero21 is offline   Reply With Quote

Reply


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 Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


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


All times are GMT -4. The time now is 10:51.