code_saturne User's Forum

Hi, is it possible set as inlet velocity a power law describing a wind speed profile?
Thanks in advance,
Regards
 

Hello,
Yes it is possible, either using the GUI with mathematical expressions (MEI), or using the usclim.f90 user subroutine.
Regards,

Thankd Yvan,

I've seen the option in the boundary conditions section.

Them I impose this law:

unorm = 0.5 * (z/10)^0.34

but I gete this error during the simulation:

SIGFPE signal (floating point exception) intercepted!
Call stack:
   1: 0x7fc6b9759bb9 <+0x5bb9>                        (libm.so.6)
   2: 0x7fc6baf74111 <+0x3111>                        (libmei.so.0)
   3: 0x7fc6baf740ed <+0x30ed>                        (libmei.so.0)
   4: 0x7fc6bb8e9d2c <uiclim_+0x222c>                 (libsaturne.so.0)
   5: 0x7fc6bba524f0 <tridim_+0x40f0>                 (libsaturne.so.0)
   6: 0x7fc6bb89cfe5 <caltri_+0x5085>                 (libsaturne.so.0)
   7: 0x7fc6bb8781cb <cs_run+0x83b>                   (libsaturne.so.0)
   8: 0x7fc6bb8784b5 <main+0x1f5>                     (libsaturne.so.0)
   9: 0x7fc6b8fcec4d <__libc_start_main+0xfd>         (libc.so.6)
  10: 0x40a259     <>                               (cs_solver)
End of stack

Whitout imposing the velocity law st the inlet, everething seems to go fine.

I'm running a steady flow simulation, should I change to unsteady?

I attach also ma xml file.

Thank you again.

Hello,
Are you using version 2.0 or 2.1 ? what are the maximum and minimum values for z ?. There might be a bug in the mathematical expression code, but to check for this, I need to build a mesh with the same range of values for z, or better, use your mesh if it is non confidential and small enough to post.
If this is the cause of the crash, using usclim.f90 while waiting for a fix would solve your immediate problem.
Another possibility is that the error occurs not when computing the prescribed velocity vectors, but when computing related turbulence inlet values (a trace similar to the one you sent with the code compiled in debug mode would tell you the line number where the crash occurs, but if I can have your mesh, I'll take a look at it anyways).
I also assume the crash occurs on the first time step ? (I hope so, as it is easier/faster to debug).
Best regards,
  Yvan

I'm using code_saturne 2.0 and here is the dropbox link to my mesh. 
yes the crash occurs to the first step.
http://dl.dropbox.com/u/5081283/Mesh_0.med
 
In the meanwhile I'll try to use usclim.f90. But can I leave my
parameters in the xml file and just set my wind profile in the subrutine?
best regards,
Filippo
 
 

Hello,
Usually, boundary conditions defined in usclim.f90 simply supersede those defined in the XMl file, and you can mix both, but I believe the XML is still fully interpreted first, so if it causes a crash, it still will.
I would recommend removing the inlet BC from the XML file, and use only that one in usclim.f90 (edit the inlet example, remove all others). All other boundary conditions should remain in the xml file.
I'll take a look at your case.
Regards,
  Yvan

Hello again,
Running under a debugger, it seems that there is a division by zero as the code tries to divide the value defined by the formula by a zero surface. This seems to be a bug in the code, and chances are (I did not test) that you could simply work around it by using "specified coordinates" instead of "normal to the inlet" as the inlet direction, using only the GUI and not usclim.f90.
I'll try to find the real cause of the crash later, and will keep you informed.
Best regards,
  Yvan

The bug is due to an spurious global initialization of a value that is later computed locally in a loop. It should be fixed in version 2.0.4 (and does not seem to occur in 2.1), but if you need a fix now, just remove lines 1405 to 1409 from src/base/cs_gui_boundary_conditions.c in the Code_Saturne kernel (ncs) source tree, and recompile/reinstall the code.
You will still get a crash a little later, as (z/10)^0.34 requires that z >= 0 (at least using the power operator on which the MEI  ^ operator is based, and for at least one face, z is negative (z = -0.00212), so using max(z, 0) instead of z should help (I have not tested this further).
Best regards,
  Yvan

thank you for the answer,
I have the possibility to install the version 2.02, does this version have the same issue or it is fixed?
Regards,
Filippo

Hello,
No, I reproduced the bug on the latest 2.0 series version (equivalent to version 2.0.3 + patches pending for 2.0.4), so you would have the same bug with version 2.0.2 or 2.0.3. The bug will be fixed in 2.0.4, but I am not sure when we will release that patch.
Otherwise, moving to 2.1 may be interesting. That version is not tested as much as 2.0, but includes many improvements (the study structure has changed a bit, so it requires re-building your case). So choosing between 2.1 and 2.0 is a matter of priority (for us developpers, feedback on the most recent version is always interesting, but for a user, priorities can be different, depending on whether you prefer the latest version, with a better data structure, better GUI, etc., but a bit less testing, or if prefer to minimize risks).
In any case, replaceing (z/10)^0.34 with (max(z,0)/10)^0.34 inyour velocity profile and using "specified coordinates" instead of "normal to the inlet" should solve your problem with your installed version (assuming the values of z below zero were due to roundoff in the meshing stage, and that z does not go significantly below 0; otherwise, you need to improve upon this).
Best regards,
  Yvan