Spatially varying body force

Questions and remarks about code_saturne usage
Forum rules
Please read the forum usage recommendations before posting.
ywan459
Posts: 20
Joined: Thu Nov 25, 2021 9:13 pm

Spatially varying body force

Post by ywan459 »

Dear all,

I am new to code_saturne and trying to setup a "simple" model for modelling the flow pattern subject to a spatially varying body force field. This body force is induced by acoustic actuator. The data is available as a series of 3D vector in a file. The data points has NO association with the computational mesh. I guess we have to go through the user subroutine or user script to read in these data and modify the existing body force fields. Could some one please point me out how to setup such facilities? Thank you very much!

Regards
ywan459
Yvan Fournier
Posts: 4070
Joined: Mon Feb 20, 2012 3:25 pm

Re: Spatially varying body force

Post by Yvan Fournier »

Hello,

Approxiy how many points are there ?
Most important, should body forces be interpolated between points or do the points represent isolated forces ?

You can add external body forces easily in code_saturne,, so the most complex part here will be projecting those forces to the computational mesh.

Best regards,

Yvan
ywan459
Posts: 20
Joined: Thu Nov 25, 2021 9:13 pm

Re: Spatially varying body force

Post by ywan459 »

Yvan Fournier wrote: Wed Feb 16, 2022 1:02 pm Hello,

Approxiy how many points are there ?
Most important, should body forces be interpolated between points or do the points represent isolated forces ?

You can add external body forces easily in code_saturne,, so the most complex part here will be projecting those forces to the computational mesh.

Best regards,

Yvan
Hi Yvan,

Yes, the body force are interpolateable. They are not isolated forces. We can certainly interpolate the force on the computational mesh through external script.

In practice, how do we add external body forces in code_saturne? I cannot find any good practical details in the user manual. Could you please give us some hits on this?

Regards
ywan459
Yvan Fournier
Posts: 4070
Joined: Mon Feb 20, 2012 3:25 pm

Re: Spatially varying body force

Post by Yvan Fournier »

Hello,

There is an "fextt" array in certain parts of the code representing external forces, but looking in detail, it does not seem directly accessible from user-defined functions (and can contain a sum of different terms).

So actually, the best way of describing the use of external forces is to define a velocity source term in cs_user_source_terms.c.

For example : https://www.code-saturne.org/documentat ... entum.html

For the interpolation part, you will need to handle that separately... If you set of points can be converted to a mesh in the MED format, and you have installed code_saturne with MEDCoupling support enabled, the code in the following example (for field initialization, but easily transposable to source terms) illustrates how that can be used for interpolation: https://www.code-saturne.org/documentat ... er_3d.html
But this does require an upstream conversion to MEDCoupling (wich could be done with SALOME's PARAVIS ParaView+plugins tool), and generation of a mesh based on your point data (for which ParaView's / PARAVIS's "Delaunay3D" filter may be useful).

Best regards,

Yvan
ywan459
Posts: 20
Joined: Thu Nov 25, 2021 9:13 pm

Re: Spatially varying body force

Post by ywan459 »

Yvan Fournier wrote: Thu Feb 17, 2022 7:35 pm Hello,

There is an "fextt" array in certain parts of the code representing external forces, but looking in detail, it does not seem directly accessible from user-defined functions (and can contain a sum of different terms).

So actually, the best way of describing the use of external forces is to define a velocity source term in cs_user_source_terms.c.

For example : https://www.code-saturne.org/documentat ... entum.html

For the interpolation part, you will need to handle that separately... If you set of points can be converted to a mesh in the MED format, and you have installed code_saturne with MEDCoupling support enabled, the code in the following example (for field initialization, but easily transposable to source terms) illustrates how that can be used for interpolation: https://www.code-saturne.org/documentat ... er_3d.html
But this does require an upstream conversion to MEDCoupling (wich could be done with SALOME's PARAVIS ParaView+plugins tool), and generation of a mesh based on your point data (for which ParaView's / PARAVIS's "Delaunay3D" filter may be useful).

Best regards,

Yvan
Hi Yvan,

I really appreciate your advise here! Thank you very much! Regarding to the velocity source term, I do not really understand how we manipulate this. The reason we add localized forces is asking solver to solve the velocity at the local region. I guess the body force field is just an independent source terms which should add to the standard NS equation. I have looked the example source term c file and found the only scaler field was added on each cell. I might be very wrong, but how do I add the vector field (having 3 components x, y and z). Sorry I am lacking knowledge of FVM implementation. I am confused here.


Regards
Yikun Wang
Yvan Fournier
Posts: 4070
Joined: Mon Feb 20, 2012 3:25 pm

Re: Spatially varying body force

Post by Yvan Fournier »

Hello,

The first example for the source terms as actually for the velocity, with a term in direction x.

And yes, the source term is an independent source term. For cases where the velocity is known at certain points, and what you want is a form of "nudging", what you need may be closer to the atmospheric imbrication module (using Cressman interpolation) for which there may be minimal examples in the code. I am not familiar with the usage of this module, so I can point to its existence if you want to search for that in the code to see if that can be useful to you, but I can't be of much more use in this case.

Best regards,

Yvan
ywan459
Posts: 20
Joined: Thu Nov 25, 2021 9:13 pm

Re: Spatially varying body force

Post by ywan459 »

Hi Yvan,

Thanks for all your helps!

I have made some progress. I have sucessfully fitted the body force data as a vector field into the med mesh file through the medcoupling interface, and my external script for data interpolation at the nodes of mesh. I found Salome's documentation is super good that I can easily follow.

Code_satune provides good examples and docs on the user-defined functions. After reading the examples c files and reference c files (came with the code_saturne installation), I have coded my source term as following which takes the acceleration field into the calculations.

Code: Select all

void
cs_user_source_terms(cs_domain_t  *domain,
                     int           f_id,
                     cs_real_t    *st_exp,
                     cs_real_t    *st_imp)
{
  /* field structure */
  const cs_field_t  *f = cs_field_by_id(f_id);

  /* mesh quantities */
  const cs_lnum_t  n_cells = domain->mesh->n_cells;
  const cs_real_t  *cell_f_vol = domain->mesh_quantities->cell_vol;
  
  /* Density */
  const cs_real_t  *cpro_rom = CS_F_(rho)->val; 
  
  if (strcmp(f->name, "VectorFieldOnNode") == 0) {
    
    cs_real_3_t    *_st_exp = (cs_real_3_t *)st_exp;
    /* apply source terms to all cells */
    for (cs_lnum_t i = 0; i < n_cells; i++) {
      _st_exp[i][0] = cell_f_vol[i] * f->vals[0][0]*cpro_rom[i];
      _st_exp[i][1] = cell_f_vol[i] * f->vals[0][1]*cpro_rom[i];
      _st_exp[i][2] = cell_f_vol[i] * f->vals[0][2]*cpro_rom[i];
    }
  }

}
The compilation is successful, and I can run the function on a toy model - a water cube surrounded by walls (6 faces are defined as walls) and subject the body force load pointing to the positive Z direction. I have compared my results with the benchmark results (defining the body force through GUI). I found my results is quite wrong. I guess something in my code must be wrong, but I cannot really figure it/them out. The major logical issue is that if

Code: Select all

 f-> vals
get the values of field, but I do not see any cell index associated with the variable. The major challenge to me are 1) whether above approach of taking the field value is all right; 2) It seems the formula of explicit term = acceleration*density*volume does not really add anything to the NS equation. I have tried to hard code the field variable as constant, and the results has no difference. Do know how what is wrong of the code. Really need professional advise here.

Thanks in advance for all your kindly helps!

Regards
ywan459
Yvan Fournier
Posts: 4070
Joined: Mon Feb 20, 2012 3:25 pm

Re: Spatially varying body force

Post by Yvan Fournier »

Hello,

Actually, you have no "solved" field in code_saturne named "VectorFieldOnNode", so the test in your function means the source term never activates.

For velocity, use
(strcmp(f->name, "velocity") == 0)
or
(f == CS_F_(vel))
which is equivalent (the second option uses predefined pointers for the main fields).

With this, you should be able to get a working source term.

The next step will be to combine this with the filtered data. I recommend testing your data input and filtering in the cs_user_postprocess section so as to be able to visualize what you are adding, and then once everything seems good, move it to source terms.

Best regards,

Yvan
ywan459
Posts: 20
Joined: Thu Nov 25, 2021 9:13 pm

Re: Spatially varying body force

Post by ywan459 »

Yes Yvan,

You are right. I changed my code to following:

Code: Select all


  if (strcmp(f->name, "velocity") == 0) {  
    cs_real_3_t    *_st_exp = (cs_real_3_t *)st_exp;
    /* apply source terms to all cells */
    for (cs_lnum_t i = 0; i < n_cells; i++) {
      _st_exp[i][0] = cell_f_vol[i] * 0*cpro_rom[i];
      _st_exp[i][1] = cell_f_vol[i] * 0*cpro_rom[i];
      _st_exp[i][2] = cell_f_vol[i] * 9.8*cpro_rom[i];
    }
  }

The results seems looking fine now (although the quantity of the total pressure does not really match to the benchmark). I misunderstood the mechanism of how to adding terms to the NS equation in code_saturne.

Reading of the spatially varying body force from "field" was not successful. Following is my code to read the data from field and add to the source term.

Code: Select all

  const cs_real_3_t *cvar_svbf = cs_field_by_name("VectorFieldOnNode")->val;
  
  /* if (strcmp(f->name, "VectorFieldOnNode") == 0) { */
  if (strcmp(f->name, "velocity") == 0) {  
    cs_real_3_t    *_st_exp = (cs_real_3_t *)st_exp;
    /* apply source terms to all cells */
    for (cs_lnum_t i = 0; i < n_cells; i++) {
      _st_exp[i][0] = cell_f_vol[i] * cvar_svbf[i][0]*cpro_rom[i];
      _st_exp[i][1] = cell_f_vol[i] * cvar_svbf[i][1]*cpro_rom[i];
      _st_exp[i][2] = cell_f_vol[i] * cvar_svbf[i][2]*cpro_rom[i];
     /*
      _st_exp[i][0] = cell_f_vol[i] * 0*cpro_rom[i];
      _st_exp[i][1] = cell_f_vol[i] * 0*cpro_rom[i];
      _st_exp[i][2] = cell_f_vol[i] * 9.8*cpro_rom[i];
      */
    }
  }

It gives the error message during the solution.

Code: Select all


/home/App/code_saturne/code_saturne-7.0.2/src/base/cs_field.c:2345: Fatal error.

Field "VectorFieldOnNode" is not defined.


Call stack:
   1: 0x7fc61b010e7f <cs_field_by_name+0x4f>          (libsaturne-7.0.so)
   2: 0x561c032801b9 <cs_user_source_terms+0x49>      (cs_solver)
   3: 0x7fc61b0fefa1 <predvv_+0x11d9>                 (libsaturne-7.0.so)
   4: 0x7fc61b0ed9f0 <navstv_+0x1214>                 (libsaturne-7.0.so)
   5: 0x7fc61b11d178 <tridim_+0x3624>                 (libsaturne-7.0.so)
   6: 0x7fc61af8cc94 <caltri_+0x1c2f>                 (libsaturne-7.0.so)
   7: 0x7fc61be58b72 <main+0x712>                     (libcs_solver-7.0.so)
   8: 0x7fc61acf60b3 <__libc_start_main+0xf3>         (libc.so.6)
   9: 0x561c032800ae <_start+0x2e>                    (cs_solver)
End of stack

It seems the "field" is only added to the mesh, but not recognized in the computational solver. I did verified the "field" in the paraview, as shown in the attached filed. I guess the solver did not picked up the field on mesh. I understand we need to define a "field" for computation in user defined parameter routine. I could not find relevant information on the manual about how to add mesh field into the computational field. The code below is how I add the mesh field through Python

Code: Select all


from MEDCouplingRemapper import *
import medcoupling as mc
import MEDLoader

import numpy as np
import csv

Path = ".//01_Data//Geometry"
file = Path + "//mesh.med"

#mesh=MEDLoader.ReadUMeshFromFile(file,"Mesh_Cube",0)
#Numofnodes = mesh.getNumberOfNodes()
#idx = list(range(0,Numofnodes,1))

mesh=MEDLoader.ReadUMeshFromFile(file,"Mesh_Cube",0)
Numofnodes = mesh.getNumberOfNodes()

#vec = np.random.randint(100, size=(Numofnodes, 3)).ravel().tolist()

vec = np.array([np.ones((1,Numofnodes))*0,
                np.ones((1,Numofnodes))*0,
                np.ones((1,Numofnodes))*9.8]).transpose().ravel().tolist()

print(vec)

fieldOnNodes=MEDCouplingFieldDouble.New(ON_NODES,CONST_ON_TIME_INTERVAL)

fieldOnNodes.setTimeUnit("s") # Time unit is seconds.
fieldOnNodes.setStartTime(0,2,-1)
fieldOnNodes.setEndTime(60000,2,-1)

fieldOnNodes.setMesh(mesh)
fieldOnNodes.setName("VectorFieldOnNode")

array=DataArrayDouble()
array.alloc(fieldOnNodes.getMesh().getNumberOfNodes(),3)

array.setValues(vec,Numofnodes,3)

fieldOnNodes.setArray(array)

MEDLoader.WriteField(Path+"//Mesh_Cube_wField.med",fieldOnNodes,True)

My question is: How to setup the field in the solver, or transferring/pick up the data from mesh field to the solver.
Attachments
Capture.PNG
ywan459
Posts: 20
Joined: Thu Nov 25, 2021 9:13 pm

Re: Spatially varying body force

Post by ywan459 »

Hi Yvan,

I think I have figured out a general workflow, the following codes are my trial setup. I would use simple example to try first, before move to the complex case. I will use 1 scalar field to define the body force universally and read out the body force in the source term.

1. define the a field in cs_user_parameters.c

Code: Select all

void
cs_user_parameters(cs_domain_t *domain)
{   /*
    cs_parameters_add_property("VectorFieldonNode",
                             3,
                             CS_MESH_LOCATION_VERTICES); */
    
    cs_parameters_add_variable("VectorFieldonNode", 1);
   
}
2. Then initialize the field variable in cs_user_initialization.c

Code: Select all

void
cs_user_initialization(cs_domain_t     *domain)
{
  const cs_mesh_t *m = domain->mesh;

  /* If this is restarted computation, do not reinitialize values */
  if (domain->time_step->nt_prev > 0)
    return;

  /* cs_property_t  *f = cs_property_by_name("VectorFieldonNode"); */
  cs_field_t *f = cs_field_by_name("VectorFieldonNode");

  if (f != NULL) {
    for (cs_lnum_t cell_id = 0; cell_id < m->n_cells; cell_id++)
      f->val[cell_id] = 9.8;
  }
}
3. Read the field variable in the cs_user_source_terms.c, in this case, we just read it out and do nothing

Code: Select all

void
cs_user_source_terms(cs_domain_t  *domain,
                     int           f_id,
                     cs_real_t    *st_exp,
                     cs_real_t    *st_imp)
{
  /* field structure */
  const cs_field_t  *f = cs_field_by_id(f_id);

  /* mesh quantities */
  const cs_lnum_t  n_cells = domain->mesh->n_cells;
  const cs_real_t  *cell_f_vol = domain->mesh_quantities->cell_vol;
  
  /* Density */
  const cs_real_t  *cpro_rom = CS_F_(rho)->val; 
  
  const cs_real_t ro0  = cs_glob_fluid_properties->ro0;
  
  /* const cs_real_t *cvar_pr = cs_field_by_name("pressure")->val; */
  /* const cs_real_3_t *cvar_svbf = cs_field_by_name("VectorFieldOnNode")->val; */
   /*cs_property_t  *pty = cs_property_by_name("VectorFieldonNode"); */
   
  const cs_real_t *cvar_svbf = cs_field_by_name("VectorFieldonNode")->val;
  
  if (strcmp(f->name, "velocity") == 0) {
    cs_real_3_t    *_st_exp = (cs_real_3_t *)st_exp;
    /* apply source terms to all cells */
    for (cs_lnum_t i = 0; i < n_cells; i++) {
      /* _st_exp[i][0] = cell_f_vol[i] * cvar_svbf[i][0]*cpro_rom[i];
      _st_exp[i][1] = cell_f_vol[i] * cvar_svbf[i][1]*cpro_rom[i];
      _st_exp[i][2] = cell_f_vol[i] * cvar_svbf[i][2]*cpro_rom[i]; */
     
      _st_exp[i][0] = cell_f_vol[i] * 0*cpro_rom[i];
      _st_exp[i][1] = cell_f_vol[i] * 0*cpro_rom[i];
      _st_exp[i][2] = cell_f_vol[i] * 9.8*cpro_rom[i];
      
    }
  }
}
But the compilation was not successful. The error message looks strange to me:

Code: Select all

/home/App/code_saturne/code_saturne-7.0.2/src/base/cs_boundary_conditions.c:699: Fatal error.


Some boundary condition definitions are incomplete or incorrect.

  For details, read the end of the calculation log,
  or visualize the error postprocessing output.


Call stack:
   1: 0x7f40880afc74 <cs_boundary_conditions_error+0xa4> (libsaturne-7.0.so)
   2: 0x7f40887fb8c5 <__cs_c_bindings_MOD_boundary_conditions_error+0xe> (libsaturne-7.0.so)
   3: 0x7f40881f10ad <vericl_+0x47fe>                 (libsaturne-7.0.so)
   4: 0x7f40880701c3 <condli_+0xf41>                  (libsaturne-7.0.so)
   5: 0x7f40881df5b8 <tridim_+0x2a64>                 (libsaturne-7.0.so)
   6: 0x7f408804fc94 <caltri_+0x1c2f>                 (libsaturne-7.0.so)
   7: 0x7f4088f1bb72 <main+0x712>                     (libcs_solver-7.0.so)
   8: 0x7f4087db90b3 <__libc_start_main+0xf3>         (libc.so.6)
   9: 0x56236ae440ce <_start+0x2e>                    (cs_solver)
End of stack
How this field could be associated with boundary conditions? I would expect no solution conducted on this field. It is mainly passing down the variables to the source term definition.

Thanks for your helps!

Regards
ywan459
Post Reply