Scaling parameters definition for electric model (cs_user_electric_scaling.cpp)

Introduction

C user function for scaling parameters definition for electric model.

Scaling parameters definition for electric model

 
  cs_lnum_t  ncel   = mesh->n_cells;
  cs_lnum_t  ncelet = mesh->n_cells_with_ghosts;
  cs_real_3_t *xyzcen =  mesh_quantities->cell_cen;
  cs_real_t *volume = mesh_quantities->cell_vol;
  cs_lnum_t  nfac   = mesh->n_i_faces;
  const cs_real_t *surfac = mesh_quantities->b_face_surf;
  const cs_nreal_3_t *u_normal = mesh_quantities->b_face_u_normal;
  const cs_real_3_t *cdgfac = mesh_quantities->i_face_cog;
 
  cs_elec_option_t *elec_opt = cs_get_glob_elec_option();
  const int kivisl = cs_field_key_id("diffusivity_id");
 
  int ielarc = cs_glob_physical_model_flag[CS_ELECTRIC_ARCS];
 
  /* example of a restrike arc */
  if (ielarc >= 1) {
    if (cs_glob_time_step->nt_cur <= 200)
      elec_opt->couimp = 200.;
    else if (cs_glob_time_step->nt_cur > 200 &&
             cs_glob_time_step->nt_cur <= 400)
      elec_opt->couimp = 200. +  2. * (cs_glob_time_step->nt_cur - 200);
    else
      elec_opt->couimp = 600.;
  }
 
  if (cs_glob_time_step->nt_cur < 400 ||
      cs_glob_time_step->nt_cur == cs_glob_time_step->nt_prev + 1)
    elec_opt->irestrike = 0;
 
  double econs = 1.5e5;
  double coepot = 0.;
  double coepoa = 1.;
 
  if (cs_glob_elec_option->irestrike) {
    double amex = 1.e30;
    double aiex = -1.e30;
    double emax = 0.;
    cs_array<cs_real_t> w1(ncelet);
    int diff_id = CS_FI_(curre, 0)->get_key_int(kivisl);
    auto cvar_prop = cs_field(diff_id)->get_val_s();
 
    auto cvar_curre0 = CS_FI_(curre, 0)->get_val_s();
    auto cvar_curre1 = CS_FI_(curre, 1)->get_val_s();
    auto cvar_curre2 = CS_FI_(curre, 2)->get_val_s();
 
    auto cvar_h = CS_F_(h)->get_val_s();
    auto cvar_rho = CS_F_(rho)->get_val_s();
 
    for (int iel = 0; iel < ncel; iel++) {
      double xelec = cvar_curre0[iel] / cvar_prop[iel];
      double yelec = cvar_curre1[iel] / cvar_prop[iel];
      double zelec = cvar_curre2[iel] / cvar_prop[iel];
 
      w1[iel] = pow(xelec * xelec + yelec * yelec + zelec * zelec, 0.5);
      amex = cs::min(amex, w1[iel]);
      aiex = cs::max(amex, w1[iel]);
    }
    cs::parall::min(amex, aiex);
 
    bft_printf("min and max for E : %14.5E %15.4E\n", amex, aiex);
 
    if (aiex > econs) {
      elec_opt->irestrike  = 1;
      elec_opt->ntdcla = cs_glob_time_step->nt_cur;
 
      /* initialize restrike point coordinates */
      elec_opt->restrike_point[0] = 1.e-8;
      elec_opt->restrike_point[1] = 1.e-8;
      elec_opt->restrike_point[2] = 1.e-8;
      double diff = 0.;
      double xyzmax[3] = {-1.e10, -1.e10, -1.e10};
 
      for (int iel = 0; iel < ncel; iel++) {
        diff = aiex - w1[iel];
 
        if (diff < 1.e-6) {
          emax = w1[iel];
          xyzmax[0] = xyzcen[iel][0];
          xyzmax[1] = xyzcen[iel][1];
          xyzmax[2] = xyzcen[iel][2];
        }
      }
      /* we can only have a single restrike point */
      cs_parall_max_loc_vals(3, &emax, xyzmax);
 
      elec_opt->restrike_point[0] = xyzmax[0];
      elec_opt->restrike_point[1] = xyzmax[1];
      elec_opt->restrike_point[2] = xyzmax[2];
 
      bft_printf("restrike point : %14.5E %14.5E %14.5E\n",
                 elec_opt->restrike_point[0],
                 elec_opt->restrike_point[1],
                 elec_opt->restrike_point[2]);
    }
 
    if (cs_glob_time_step->nt_cur <= elec_opt->ntdcla + 30) {
      double z1 = elec_opt->restrike_point[0] - 3.e-4;
      double z2 = elec_opt->restrike_point[0] + 3.e-4;
      if (z1 < 0.)
        z1 = 0.;
      if (z2 > 2.e-2)
        z2 = 2.e-2;
      for (int iel = 0; iel < ncel; iel++) {
        if (xyzcen[iel][2] > z1 && xyzcen[iel][2] < z2) {
          double rayo =   elec_opt->restrike_point[0] * xyzcen[iel][0]
                        - elec_opt->restrike_point[1] * xyzcen[iel][1];
          double denom
            = pow(  elec_opt->restrike_point[0] * elec_opt->restrike_point[0]
                  + elec_opt->restrike_point[1] * elec_opt->restrike_point[1], 0.5);
          rayo /= denom;
          rayo += (xyzcen[iel][2] - elec_opt->restrike_point[2]
                 * xyzcen[iel][2] - elec_opt->restrike_point[2]);
          rayo = pow(rayo, 0.5);
 
          double posi = elec_opt->restrike_point[0] * xyzcen[iel][0];
 
          if (rayo < 5.e-4 && posi <= 0.)
            cvar_h[iel] = 8.e7;
        }
      }
    }
    else {
      elec_opt->irestrike = 0;
    }
    double somje = 0.;
 
    auto cvar_joulp = CS_F_(joulp)->get_val_s();
    for (int iel = 0; iel < ncel; iel++) {
      somje += cvar_joulp[iel] * volume[iel];
    }
 
    cs::parall::sum(somje);
 
    if (fabs(somje) > 1.-20)
      coepot = cs_glob_elec_option->couimp * cs_glob_elec_option->pot_diff
              / cs::max(somje, cs_math_epzero);
 
    bft_printf("imposed current %14.5E, Dpot %14.5E, Somje %14.5E\n",
               cs_glob_elec_option->couimp,
               cs_glob_elec_option->pot_diff,
               somje);
 
    double elcou = 0.;
    for (int ifac = 0; ifac < nfac; ifac++) {
      if (fabs(u_normal[ifac][0]) < 1.e-2 && fabs(u_normal[ifac][1]) < 1.e-2 &&
          cdgfac[ifac][2] > 0.05e-2 &&  cdgfac[ifac][2] < 0.08e-2) {
        int iel = mesh->i_face_cells[ifac][0];
        elcou += cvar_curre2[iel] * u_normal[ifac][2] * surfac[ifac];
      }
    }
 
    cs::parall::sum(elcou);
 
    if (fabs(elcou) > 1.e-6)
      elcou = fabs(elcou);
    else
      elcou = 0.;
 
    if (fabs(elcou) > 1.e20)
      coepoa = cs_glob_elec_option->couimp / elcou;
 
    coepot = coepoa;
 
    double dtj = 1.e15;
    double dtjm = dtj;
    double delhsh = 0.;
    double cdtj = 20.;
 
    for (int iel = 0; iel < ncel; iel++) {
      if (fabs(cvar_rho[iel]) > 1.e-20)
        delhsh = cvar_joulp[iel] * dt[iel]
               / cvar_rho[iel];
 
      if (fabs(delhsh) > 1.e-20)
        dtjm = cvar_h[iel] / delhsh;
      else
        dtjm = dtj;
      dtjm = fabs(dtjm);
      dtj = cs::min(dtj, dtjm);
    }
    cs::parall::min(dtj);
 
    double cpmx = pow(cdtj * dtj, 0.5);
    coepot = cpmx;
 
    if (cs_glob_time_step->nt_cur > 3) {
      if (coepoa > 1.05)
        coepot = cpmx;
      else
        coepot = coepoa;
    }
 
    bft_printf(" Cpmx          = %14.5E\n", cpmx);
    bft_printf(" COEPOA        = %14.5E\n", coepoa);
    bft_printf(" COEPOT        = %14.5E\n", coepot);
    bft_printf(" Dpot rescaled = %14.5E\n",
               cs_glob_elec_option->pot_diff * coepot);
 
    /* scaling electric fields */
    elec_opt->pot_diff *= coepot;
 
    /* electric potential (for post treatment) */
    auto cvar_potr = CS_F_(potr)->get_val_s();
    for (int iel = 0; iel < ncel; iel++)
      cvar_potr[iel] *= coepot;
 
    /* current density */
    if (ielarc > 0) {
      for (int iel = 0; iel < 3 ; iel++) {
        cvar_curre0[iel] *= coepot;
        cvar_curre1[iel] *= coepot;
        cvar_curre2[iel] *= coepot;
      }
    }
 
    /* joule effect */
    for (int iel = 0; iel < 3 ; iel++)
      cvar_joulp[iel] *= coepot * coepot;
  }