Initialization from a 3D post-processing output using MEDCoupling

Use Case

When CFD results from another tool are available, they may be used to initialize a code_saturne computation using the MEDCoupling remapper tools. For more information on MEDCoupling features, see https://docs.salome-platform.org/latest/dev/MEDCoupling/tutorial/index.html#english-toc

This requires that:

• code_saturne is installed with MEDCoupling
• MEDCoupling is installed with MED file support (which is the default configuration).
• The provided data is in MED format or at least in a format supported by ParaView.

If the previous results are provided in another format than MED, they can be converted using the SALOME platform's PARAVIS module, which is an extension to ParaView including several plugins, especially a MED reader and writer. All that is required is to read the results in PARAVIS, then save them to MED format them using the "Save Data" option from ParaView's "File" menu. This can also be done using ParaView scripts if preferred.

Example code

The following example shows how to read fields from a MED file using MEDCoupling remapper features.

The field_names array here contains the list of fields we plan to read and their matching names in the MED file.

The file_name should reference a file relative to the run directory, using either a proper relative path or an absolute path (as usual, if the file is placed in a case's DATA directory, it will automatically be copied to the execution directory, but this leads to additional copies of possibly large data).

In this example, we also force the remapper option to PointLocator (see INTERP_KERNEL::PointLocator in MEDCoupling documentation). By default, code_saturne uses INTERP_KERNEL::Triangulation, which allows for better conservation of quantities, but is slower.

  /* Apply only on computation start, not on restarts */
  if (domain->time_step->nt_prev > 0)
    return;

#if defined(HAVE_MEDCOUPLING_LOADER)

  double t1 = cs_timer_wtime();

  const cs_lnum_t n_cells = cs_glob_mesh->n_cells;

  bft_printf("In cs_user_initialization HAVE_MEDCOUPLING_LOADER\n");

  /* Number of fields to interpolate from the MED file */
  const int  nremapper_fields = 2;

  /* Names of the nremapper_fields fields to read */
  const char  *field_names[] = {"p", "U"};

  /* We request a remapper with a given name. If it does not exist,
   * the function returns a NULL pointer. */
  cs_medcoupling_remapper_t *r
    = cs_medcoupling_remapper_by_name_try("init");

  if (r == NULL) {

    /* Space dimension of the elements (2 for faces, 3 for cells) */
    int elts_dim = 3;

    /* Indexes needed to read the time step from the
     * file (0, 0 if only one exists) */
    int it0 = 0, it1 = 0;

    /* Path to file */
    const char file_name[] = "../../../MESH/channel395_OF_4.1_5000.med";

    /* The remapper is created. We retrieve its id from the function.
       * The function inputs are:
       * 1) Name of the remapper
       * 2) dimension of the mesh elements
       * 3) selection criteria
       * 4) path to the med file
       * 5) number of fields to interpolate
       * 6) names of the fields to interpolate
       * 7 + 8) time iteration index and order */
    int r_id = cs_medcoupling_remapper_initialize("init",
                                                  elts_dim,
                                                  "all[]",
                                                  file_name,
                                                  nremapper_fields,
                                                  field_names,
                                                  it0,
                                                  it1);

    /* Retrieve the pointer */
    r = cs_medcoupling_remapper_by_id(r_id);

    /* Set options */
    cs_medcoupling_remapper_set_options(r,
                                        "IntersectionType",
                                        "PointLocator");

    /* We create the interpolation matrix => Here it is only called once
     * since the mesh is not moving */
    cs_medcoupling_remapper_setup(r);

  }

  /* We retrieve arrays containing the interpolated values.
   * Inputs are:
   * 1) remapper
   * 2) id of the field to interpolate
   * 3) a default value (if no intersection is obtained) */

  /* We loop in the fields matching field_names defined above */

  for (int f_id = 0; f_id < nremapper_fields; f_id++) {

    cs_real_t *m_vals = cs_medcoupling_remapper_copy_values(r, f_id, 1.);

    /* For pressure */

    if (f_id == 0) { /* or (strcmp(field_names[f_id], "p") == 0) */

      cs_real_t *cvar_p = CS_F_(p)->val;

      for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
        cvar_p[c_id] = m_vals[c_id];
      }

    }

    /* For velocity */

    else if (f_id == 1) { /* or (strcmp(field_names[f_id], "U") == 0) */

      cs_real_3_t *cvar_vel = (cs_real_3_t *)CS_F_(vel)->val;

      for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
        cvar_vel[c_id][0] = m_vals[c_id*3];
        cvar_vel[c_id][1] = m_vals[c_id*3 + 1];
        cvar_vel[c_id][2] = m_vals[c_id*3 + 2];
      }
    }

    BFT_FREE(m_vals);

  } /* End of loop on fields */

  double t2 = cs_timer_wtime();

  bft_printf("Time in MEDCoupling: %f seconds \n", t2 - t1);
  bft_printf_flush();

#endif

Note that we do not need to specify the field dimensions here, but they should match, so mapping a vector field to a scalar will produce a crash.