Lagrangian simulations - Additional variables

[Pablo]

Hello everybody:

This is my first post both in this board and in Code_Saturne community so, first of all I think a brief introduction should be mandatory:

My name is Pablo and I am performing my PhD with Code_Saturne since a (relatively) recent time. Because of my academic/professional past I have some experience with other commercial softwares but since I discovered Code_Saturne (a few years ago) I have boosted its use and knowledge due to it is a unfairly unknown Open Source tool out of research circles.

The point is my PhD concerns Lagrangian simulations and so I am developing preliminary simulations which involves additional variables which are calculated as function of other variables and post-processed. By now, the strategy I am following is based on:

- Frozen-field preliminar simulation to define the velocity/thermal field.
- Lagrangian simulation based on the aforementioned simulation defined by GUI/subroutines:
- Subroutine "uslag1":Allows to define the additional variable "nvls=1", as well as other simulation parameters
- Subroutine "uslaed": Integration of the additional variable SDE
- Subroutine "cs_user_extra_operations": Allows to extract, calculate and save the additional variable in the "ettp" array.

The question is, although the calculation and save process in the "ettp" array is correctly made (because the values stored in the "ettp(npt, jvls(1))" has been checked with a "write(nfectra,...)" statement), I can't find a way to check the results in the postprocesor (Paraview) because I haven't been able to find in the documentation a way to "record" the additional variable on the particle trajectories as for example temperature would be recorded via the activation of the "iviste=1" in the "uslag1" subroutine.

Could somebody help me to visualize the additional variable results on the particle trajectories?

Thanks in advance.



PS. I'm actually using both Code_Saturne 3.2 / 3.2.2 versions

[Yvan Fournier]

Hello,

The use of extra Lagrangian variables has been broken in Code_Saturne 3.0 (due to the rewrite of the trajectography for parallelism, which was not quite complete).

I "repaired" this with some improvements to the Lagrangian module last spring, meaning ettp(npt, jvls(1)) should work again, but only starting from version 3.3. (It has no chance of working in versions 3.0 to 3.2). As it is a rarely used feature at EDF, I did not test this feature in 3.3 (it works in theory), so if you switch to 3.3, please let me know if it works, and if it does not, I'll be happy to have a test case to debug it.

For postprocessing, you need to use the "cs_post_write_particle_values" function, for which there is currently no Fortran wrapper. I'll try to prepare one for 4.0 (which, as a "long term support"/fully validated version, will be branched for Beta sometimes in October, which is arriving fast). For the meantime, although I dont' usually recommend modifying non-user functions, the easiest is probably for you to copy the base/src/cs_post_default.c file from the sources to your case's SRC directory, and adapt the _write_particle_vars function in there. CS_LAGR_N_ATTRIBUTES is the number of possible non-user variables. Your user variables should come just after that, so adding the following after the loop (to avoid the test on input->particle_attr[attr_id], which is only valid for non-user values) should work:

Code: Select all

cs_post_write_particle_values(mesh_id,
                              CS_LAGR_N_ATTRIBUTES + user_attr_id, /* user_attr_id starts at 0*/
                              "user_var_1", /* your choice here */
                              0, /* component_id */
                              ts);

Regards,

Yvan

[Pablo]

Hello Yvan:

First of all, thanks you for the quick answer to this doubt.

I guess from your answer that the proposed solution (the one which implies the addition of several group of lines to the cs_post_default.c file and its copy to the simulation SRC directory) could only work in Saturne 3.3, which implies I have to update Code_Saturne. No problem about that, I will update the program although I will test it anyway considering it is not going to work (in Saturne 3.2.2, which is the version I am testing right now).

Anyway, there is a few comments I would like to be checked by you (if possible) considering these tests has been perfomed on the Code_Saturne 3.2.2:

a) You said you repaired the Lagrangian module which allows the ettp(npt,jvls(1))to work again. Well, my tests indicates this command is working right, as the following lines, in the cs_user_extra_operations subroutine:

!Variables initialization
usp1=0.d0
usp2=0.d0

do npt=1,nbpart

usp1=ettp(npt,jvp)
usp2=5 * usp1
ettp(npt,jvls(1))=usp2
usp3=ettp(npt,jvls(1))

write(nfecra,90108) usp1, usp2, usp3
90108 format(&
'------------------------',/,&
'INFORMATION: ',/,&
'Variable no1 = ',f10.3 ,/,&
'Variable no2 = ',f10.3 ,/,&
'Variable no3 = ',f10.3 ,/,&
'-',/)
enddo

...seems to work so that a Lagrangian property is extracted (usp1=ettp(npt,jvp)), then calculated (usp2=5 * usp1), stored in the Lagrangian variables array (ettp(npt,jvls(1))=usp2), extracted again (usp3=ettp(npt,jvls(1))) and shown in the "output window" while calculating with the write command. This indicates the ettp(npt,jvls(1)) works correctly despite it is being ran on Code_Saturne 3.2.2, doesn't it?
(I would like to check if the input/output commands on the Lagrangian variables arrays are being used correctly)

b) Considering your comments, I have re-scripted the _write_particle_varsfunction in the cs_post_default.c file (already copied on the CASE/SRC folder) as it follows:


_write_particle_vars(cs_post_default_input_t *input,
int mesh_id,
const cs_time_step_t *ts)
{
cs_lagr_attribute_t attr_id;

char var_name[64];
int component_id;
char var_name_component[64];

for (attr_id = 0; attr_id < CS_LAGR_N_ATTRIBUTES; attr_id++){

if (input->particle_attr[attr_id]) {

/* build name */

int i;
int l = snprintf(var_name,
63,
"particle_%s",
cs_lagr_attribute_name[attr_id] + strlen("cs_lagr_"));
var_name[63] = '\0';
for (i = 0; i < l; i++)
var_name = tolower(var_name);

/* Output values */

if (input->particle_multicomponent_export[attr_id] == -1)
cs_post_write_particle_values(mesh_id,
attr_id,
var_name,
input->particle_multicomponent_export[attr_id],
ts);
else {
/* Create one output per component */
for (component_id = 0; component_id < input->particle_multicomponent_export[attr_id]; component_id++) {
snprintf(var_name_component,
63,
"%s_layer_%2.2i",
var_name,
component_id+1);
var_name_component[63] = '\0';
cs_post_write_particle_values(mesh_id,
attr_id,
var_name_component,
component_id,
ts);
}
}
}

}
/*------------------------------------------------------------------------------------*/
/*--- LINE MODIFICATION TO ALLOW USER VARIABLES RESULTS ON PARTICLES TRAJECTORIES ---*/

cs_post_write_particle_values(mesh_id,
CS_LAGR_N_ATTRIBUTES + 0, /*user_attr_id starts at 0*/
"user_var_1", /*your choice here*/
1, /* component_id*/
ts);

/*------------------------------------------------------------------------------------*/

}

...which I think is the correct position as you indicated (adding a couple of extra brackets to ensure the new lines are OUT of the for (attr_id = 0; attr_id < CS_LAGR_N_ATTRIBUTES; attr_id++) loop). Is it your suggested locationfor the extra lines? I have a couple of doubts:

- user_attr_id starts at 0, so, if I have defined only one extra variable at lagrangian array position jvls(1), the extra variable will have a user_attr_id=0 (or user_attr_id=1?)
- Is the component_id related with the particle "class"?, is it required to be indicated as "1" as only one particle class is defined, or is it valid if it is indicated as the variable istelf"component_id" instead of value "1"?

Thanks very much for your attention, I will keep you informed of the Saturne 3.3 results.


Regards,


Pablo.

[Pablo]

Hello Yvan:

I have installed the Code_Saturne 3.3 from scratch in an Ubuntu 14.04, and I am not able to find the cs_post_default.c file, so I can't find where should I write the extra lines in order to be able to compile it and so, test it.

May the function cs_post_write_particle_values function placed in other file?
Also, I have seen that the files distribution in the aforementioned installation (3.3) is rather different than the previous installation (3.2.2) which was made in a CAElinux distribution.


Thanks for your replies.

Pablo

EDIT: Ok, it's my fault . Saturne 3.3 was not correctly installed. Now it's correctly installed I've been able to find the "cs_post_default.c" file and so, be able to apply the suggested solution. I will inform you soon about the results

[Yvan Fournier]

Hello,

How "hard" was your test for the user values which "seem to work" ? Did you test with constant values, or something which changes with time ? The issue is that in version 3.0 to 3.2, the values placed in ettp(npt,jvls(1)) were not handled in the trajectory stage, so they are probably not updated during that stage.

Meaning that if you simply recompute the values at each time step after the particle displacement stage, you code probably in fact does work on 3.2, but if you need to transport quantities, it should not work correctly...

Regarding your postprocessing code, it seems OK at quick glance (I did not check it in detail). I may also have a better solution, using only cs_user_postprocess.c, so I may send you another example to try in a few days, but I need to think about it first.

The extra variable will indeed have usr_attr_id = 0. Also, if I remember correctly the component id may be used both for multiple components (such as x,y,z for position and velocity), or for multiple values in subarrays, such as those quantities for multiple coal particle "layers" representing a particle as a sphere with non-constant properties between the center and the outside. Here also, they start at 0.

To test you visualization code, I recommend first setting a simple value in you user variable, such as the current x or y coordinate. That way, a quick visualization will tell you if the visualized values are OK (expected) or not (unexpected).

Regards,

Yvan

[Pablo]

Hello Yvan:

Regarding your questions, I answer:

The tests has been perfomed on the velocity values of each particle, in two different ways:
- First: Read the value of "y" direction velocity for each particle, store in an auxiliary variable and show it in the "calculation screen" (which makes the calculation go slower....) The test was performed succesfully so that the values per time step were checked to be changing with time.
- Second: Read the value of "y" direction velocity for each particle, store in an auxiliary variable and modify the z-axis velocity value for each particle with its result. This test was successful too because the particles behaviour in the post-processor was altered as it was supposed to be done.

...based on these tests, it seems the "ettp" command is working correctly so that existing particle properties may be changed as a certain calculation indicates with time. Anyway, I will try to alterate a certatin particle property so that its value will be constant in order to check if its stored in time correctly, but it seems it will works.

On the other hand, I have already tested the suggested code in CS 3.3 and although there aren't errors, results are still not visualized in the postprocess (both cs_user_extra_operations.f90 and CS_3_3_cs_post_default.c are attached with this post), but several writing errors during simulation execution happens:

...if I define the attr_id so that its value is CS_LAGR_N_ATTRIBUTES + 0, an error in the simulation happens with a text like:

cs_lagr_extract.c:325: Fatal error.

Attribute 35 has a number of components equal to 1
but component 1 is requested.

...which I think it means the CS_LAGR_N_ATTRIBUTES + 0 = attr_id = 35 position is occupied, as this error happens again when CS_LAGR_N_ATTRIBUTES + 1 = attr_id = 36. BUT, it does not happens when CS_LAGR_N_ATTRIBUTES + 2 = attr_id = 37.

Does this means the attr_id = 37 is free to be filled with the jvls(1) additional property via component_id identification?
but...
...I still haven't found how to identify the component_id which is related to the jvls(1) property which is being calculated and stored in the cs_user_extra_operations.f90 as I haven't been able to find information about this of any kind ...

Sorry for inconveniences, I think this might be a very specfic error hard to be solved until the next CS release.
Thanks for your attention.

[Yvan Fournier]

Hello,

Did you try with attr_id = 0 and component_id < 0 ? (or componen_dt = -1) ?

Component ids also start at 0, and id -1 means all components (same as 0 for a single component, all components for vectors).

Regards,

Yvan

[Pablo]

Hello Yvan:

Thanks to your suggestion now the simulation is able to store and represent the values of a different result differente than the typical ivisv1, ivisv2, ivisdm,... etc., but to be honest, I think I am lost with what is being represented or why are the values so different to what it was expected to be represented.

I have changed the values of the lagrangian additional variable to a constant value (value = "81" to all the values of the additional variable jvls(1) ) in order to ease its identification when the postprocess is done, but the range of results when the postprocess is done is far from the expected results (for instance, when attr_id = 0 and componen_dt = -1 results varies from 861 to 60133 in the first time step and from 129 to 67814 in the last one)

On the other hand, answering in your post:

Hello,

Did you try with attr_id = 0 and component_id < 0 ? (or componen_dt = -1) ? Pablo: attr_id<1 means erros in the simulation stage with a SIGSEGV signal (forbidden memory area access) intercepted! message. Simulation only works with attr_id = 0

Component ids also start at 0, and id -1 means all components (same as 0 for a single component, all components for vectors).Pablo: Unexpected results happens wether componen_dt = -1 or componen_dt = 0 is used (componen_dt = 1, 2 or -2 makes fail the simulation stage), but with same post process results for a theoretical constant value of "81"

Regards,

Yvan

...so my questions now are related with the post-process treatment of the _write_particle_vars function:

- If attr_id is a variable index of a result array with a maximun value of "CS_LAGR_N_ATTRIBUTES" and considering the existent loop manages "standard calculated" variables from attr_id=0 to attr_id = CS_LAGR_N_ATTRIBUTES (prior modification in the cs_post_default.c file), what is the point on storing a attr_id < CS_LAGR_N_ATTRIBUTES variable (in which range the additional lagrangian variable is not supposed to be, and you have suggested to use attr_id=0 )?

- Is thisfunction storing the results per time-step and particle or the results are being "summed" per particle as time-steps goes by?

- I haven't understant your comment concerning component id = 0 is related to "single component" and -1 to "all components", if, as far as I know, the additional lagrangian variable has not "directional character" related to a vector being just a calculated value based on another lagrangian results (or, that is what it will be supposed to be done in a certain future).


Once again, thank you very very much for your help, and sorry for not being enough clear in the descriptions and/or technical program understanding to allow you help me

[Yvan Fournier]

Hello,

As the functionnality you are using is not tested much, it is hard to guess whether there is in error in your code or in our code... If you have a small test case representing that you could post (mesh, data, and sources, including your latest postprocessing attempt), I could try to run it and debug...

Best regards,

Yvan

[Yvan Fournier]

Hello,

Working on the Lagrangian code yesterday, I realized I made a mistake in my suggestions.

The CS_LAGR_USER attribute is already present, so in all my previous posts on this thread, you can replace CS_LAGR_N_ATTRIBUTES by CS_LAGR_USER.

This should help things work better...

Sorry for the wrong directions.

Regards,

Yvan