access fluid properties, cs_wall_function.h

[Ionut G]

Hello,

I am new to code_saturne, I want to implement a different wall function than the ones that are already available for the k-e turbulence model.

From what I understood it is simple just to replace, in cs_wall_functions.h file, the equations for a built-in wall function with the equations required by the wall function that I want to implement.

For the wall function I have, I need to access the density, molecular viscosity, wall shear stress and the pressure gradient.
I tried to use "const cs_real_t rho = CS_F_(rho)->rho" and "const cs_real_t mu = CS_F_(mu)->mu", but it does not work.
I receive the following warnings:
warning: implicit declaration of function ‘CS_F_’ [-Wimplicit-function-declaration]
warning: nested extern declaration of ‘CS_F_’ [-Wnested-externs]
I attached the compile.log file containing the warnings

1. How can I fix these warnings?
2. How can I calculate the pressure gradient at the cells? Someone recommended I should calculate it in the clptur.f90 file. I found a subroutine "gradient_s" that is supposed to calculate the gradient of a variable, but I do not know how to use it.
Is there an example or documentation on how to use it?

Thanks,
Ionut

Code: Select all

inline static void
cs_wall_functions_2scales_log(cs_real_t   l_visc,
                              /* cs_real_t   t_visc, */
                              cs_real_t   vel,
                              cs_real_t   y,
                              /* cs_real_t   kinetic_en, */
                              int        *iuntur,
                              cs_lnum_t  *nsubla,
                              cs_lnum_t  *nlogla,
                              cs_real_t  *ustar,
                              cs_real_t  *uk,
                              cs_real_t  *yplus,
                              cs_real_t  *ypup,
                              cs_real_t  *cofimp)
{
  /* const double ypluli = cs_glob_wall_functions->ypluli; */
	
	/* Local variables */
	double wss, spg, utau, up, alpha;
	int sign_wss, sign_spg; /* signs of wss and spg */
	
        const cs_real_t rho = CS_F_(rho)->rho;
	const cs_real_t mu = CS_F_(mu)->mu;
	
	/* Compute Wall Shear Stress, wss */
	wss = 0.001003 * (vel / y);
	if (wss < 0.)				{ sign_wss = -1; }
	else if (wss > 0.)			{ sign_wss = 1;   }
	else 						{ sign_wss = 0;  }

	/* Compute Streamwise Pressure Gradient, spg */
	spg = gradpr;
	if (spg < 0.)				{ sign_spg = -1; }
	else if (spg > 0.)			{ sign_spg = 1;   }
	else						{ sign_spg = 0;  }

	/* Compute the characteristic velocity, function of wss and spg */
	utau = pow(fabs(wss) / 998.2, (1. / 2.)); /* velocity due to wss */
	up = pow((0.001003 / pow(998.2, 2)) * fabs(spg), (1. / 3.)); /*velocity due to spg*/
	*uk = sqrt(pow(utau, 2) + pow(up, 2)); /*characteristic velocity, utaup*/
	
	
	alpha = pow(utau, 2.) / pow(*uk, 2.); /* separation point indicator */
	*yplus = (*uk * y) / l_visc; /*dimensionless distance, ystar*/
	
	/* Manhart wall function model, dimensionless velocity */
	*ustar = sign_wss * alpha * *yplus +
		       sign_spg * 0.5 * pow((1 - alpha), 1.5) * pow(*yplus, 2);
				 
	/* Viscous sub-layer */
	*ypup = 1.;
	*cofimp = 0.;
	*iuntur = 0;
	*nsubla += 1;
	
	/* Viscous sub-layer */
	*nlogla = 0.;
	
}

[Yvan Fournier]

Hello,

Since the relevant code is in a header (.h) file, to ensure modifications are accounted for, you must ensure the code is recompiled using your modified version, or all source files calling these functions are also copied in the SRC folder of a case. Otherwise, only the sources in your SRC folder will use you modified version. I usually recommend the first option, but since you seem to be running on Windows, you might only have the second option.

The warnings relative to vscanf and snprinf are due to system header inclusion options, so there may not be much do do about those (at least testing slight changes in system definition macros is not very practical and might not be sufficient).

The "implicit declaration of function ‘CS_F_’ [" warning means you need to include "cs_field_pointer.h" in your modified file. Also, gradpr is an array, not a function. Computing the gradient for each face is not efficient, and not even possible in parallel mode (and the cs_wall_functions_2scales_log function is called for each wall face), so the simplest solution is probably to compute the gradient once per time step (for example in cs_user_initialize and cs_user_extra_operations.c), placing it in a global array (allocated in cs_user_initialize), and referencing it from cs_wall_functions_2scales_log).

Also, I do not understand how you determine whether the pressure gradient is streamwise or upstream, as you test uses only the sign but no 3D geometry projections or trigonometry functions.

Regards,

Yvan

[Ionut G]

Thank you very much, Yvan, for support.

I am thinking that from the field velocity I can calculate a unit vector in the direction of the flow, and based on this unit vector I can determine if the pressure gradient is streamwise or not.

The streamwise pressure gradient is not calculated yet in my code.

How can I calculate the pressure gradient near the wall?
I saw in another post that you recommended that it should be used "gradient_s" instead of "field_gradient_scalar", which one is best suited to be used?
viewtopic.php?f=2&t=2256&p=12334&hilit= ... t_s#p12334

Regards,
Ionut

[Yvan Fournier]

Hello,

I recommended "gradient_s" in the case of non-resolved variables (because they lack boundary condition definitions), but as pressure is a resolved variable, using "field_gradient_scalar" is a better option here, as it will produce a result consistent with the defined pressure boundary conditions.

Regards,

Yvan

[Ionut G]

Thank you very much, Yvan.

I wrote a code in "cs_user_initialization.f90" for calculating the streamwise pressure gradient.
To access the pressure gradient I have used the "field_gradient_scalar" subroutine and seems to work fine.

When I try to calculate the magnitude of the velocity vector I receive the following error:
"Error: Unexpected STATEMENT FUNCTION statement at (1)"
How can I fix this issue?
I attached the code and the compile log.

Regards,
Ionut

[Yvan Fournier]

Hello,

You declared "mag_vel" as a single "double precision" float, but ar using it as an array. Similarly, spg does not seem declared. This explains the compilation issue.

If you need to use those from another function, declare them in cs_user_modules.f90 (do not forget to allocate/initialize them), and add "use cs_user_modules" in your function. Otherwise, simply declare them and allocate them here.

Regards,

Yvan

[Ionut G]

Thank you Yvan, the code is working now without errors.

Now I have two problems:
1.
I want to use the "spg" variable in the cs_wall_functions.h file.
How can I call the "spg" variable from cs_wall_functions.h file?

2.
I saw that if I use 1 CPU core the simulation is fine, but when I try to use more CPU cores I receive the following error:
SIGSEGV signal (forbidden memory area access) intercepted!

Code: Select all

****************************************

  Compiling user subroutines and linking

 ****************************************

 ****************************

  Preparing calculation data

 ****************************

 Parallel Code_Saturne on 3 processes.

 ***************************

  Preprocessing calculation

 ***************************

 **********************

  Starting calculation

 **********************

job aborted:

[ranks] message

[0-1] terminated

[2] application aborted

aborting MPI_COMM_WORLD (comm=0x44000000), error 1, comm rank 2

---- error analysis -----

[2] on Ionut-PC

.\cs_solver.exe aborted the job. abort code 1

---- error analysis -----

 ****************************

  Saving calculation results

 ****************************

solver script exited with status 1. 

Error running the calculation. 

Check Code_Saturne log (listing) and error* files for details. 

Error in calculation stage. 

How can I fix this error and make it run on multiple cores?

Regards,

Ionut

[Yvan Fournier]

Hello,

I do not know where/how you defined the "spg" variable, so cannot tell you how to access it.

For the parallel crash, this is probably a bug in your code, but since we do not have a stack trace on Windows, I can't be sure (debugging on Linux would be easier).

Regards,

Yvan

[Ionut G]

Hi,

The "spg" variable is defined in cs_user_initialization-spg.f90 file, I attached it.
I want to use the values of "spg" variable in the cs_wall_functions.h, in the following part of the code:
In the cs_wall_functions.h file I defined the "spg" variable as double. I attached the cs_wall_functions.h/.c too.

Code: Select all

#include "cs_field_pointer.h"

inline static void
cs_wall_functions_2scales_log(cs_real_t   l_visc,
                              cs_real_t   vel,
                              cs_real_t   y,
                              int        *iuntur,
                              cs_lnum_t  *nsubla,
                              cs_lnum_t  *nlogla,
                              cs_real_t  *ustar,
                              cs_real_t  *uk,
                              cs_real_t  *yplus,
                              cs_real_t  *ypup,
                              cs_real_t  *cofimp)
{
	/* Fluid properties */
	const cs_real_t *rho = CS_F_(rho)->val;
	const cs_real_t *mu = CS_F_(mu)->val;
	/* Local variables */
	double wss, spg, utau, up, alpha;
	int sign_wss, sign_spg; /* signs of wss and spg */
	/* Compute Wall Shear Stress, wss */
	wss = *mu * (vel / y);
	if (wss < 0.)		{	sign_wss = -1; }
	else if (wss > 0.)	{	sign_wss = 1;	 }
	else 				{	sign_wss = 0;	 }

	/* Compute Streamwise Pressure Gradient, spg */
	
	spg; /* The value imported from cs_user_initialization-spg.f90 file */
	
	if (spg < 0.)		{ sign_spg = -1; }
	else if (spg > 0.)	{ sign_spg = 1;   }
	else				{ sign_spg = 0;  }
	
	/* Compute the characteristic velocity, function of wss and spg */
	utau = pow(fabs(wss) / *rho, (1. / 2.)); /* velocity due to wss */
	up = pow((*mu / pow(*rho, 2)) * fabs(spg), (1. / 3.)); /*velocity due to spg*/
	*uk = sqrt(pow(utau, 2) + pow(up, 2)); /*characteristic velocity, utaup*/
	alpha = pow(utau, 2.) / pow(*uk, 2.); /* separation point indicator */
	*yplus = (*uk * y) / l_visc; /*dimensionless distance, ystar*/
	/* Manhart wall function model, dimensionless velocity */
	*ustar = sign_wss * alpha * *yplus + sign_spg * 0.5 * pow((1 - alpha), 1.5) * pow(*yplus, 2);		 
	/* Viscous sub-layer */
	*ypup = 1.;
	*cofimp = 0.;
	*iuntur = 0;
	*nsubla += 1;
	/* Viscous sub-layer */
	*nlogla = 0.;
}

Thank you,
Ionut

[Yvan Fournier]

Hello,

Defining an array in Fortran and accessing it from C is a bit involved, but in you case, the simplest solution is probably to define a cell-based field in cs_user_parameters.c or cs_user_parameters.f90 (see Doxgen examples).

For example, in cs_user_model, in cs_user_parameters.c, add:

Code: Select all

cs_parameters_add_property("spg",
                           1,
                           CS_MESH_LOCATION_CELLS);

Field values are accessible both from C and from Fortran, so then in cs_user_initialization.c:

Code: Select all

cs_real_t  *cs_glob_val_spg; /* defined outside functions so as to be global */

void
cs_user_initialization(void)
{
  cs_field_t *f = cs_field_by_name("spg");
  cs_glob_val_spg = f->val;
}

Note that this stage cannot be done in cs_user_parameters.c, as the field values are only allocated after the mesh is read (after cs_user_model is called).

Finally, in cs_wall_functions.h, you can add:

Code: Select all

extern cs_real_t *cs_glob_val_spg;

and use cs_global_val_spg as your spg array (you could also name it "spg" by I recommend using name prefixes so as to emulate namespaces for global variables, to avoid possible name collisions.

Regards,

Yvan