diff --git a/src/Kernel/VF/Champs/Champ_front_softanalytique.cpp b/src/Kernel/VF/Champs/Champ_front_softanalytique.cpp index 1904c040d6..59edeb2adf 100644 --- a/src/Kernel/VF/Champs/Champ_front_softanalytique.cpp +++ b/src/Kernel/VF/Champs/Champ_front_softanalytique.cpp @@ -1,5 +1,5 @@ /**************************************************************************** -* Copyright (c) 2022, CEA +* Copyright (c) 2026, CEA * All rights reserved. * * Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: @@ -16,6 +16,7 @@ #include #include #include +#include Implemente_instanciable(Champ_front_softanalytique,"Champ_front_fonc_xyz",Ch_front_var_stationnaire); // XD champ_front_fonc_xyz front_field_base champ_front_fonc_xyz 0 Boundary field which is not constant in space. @@ -93,39 +94,53 @@ int Champ_front_softanalytique::initialiser(double temps, const Champ_Inc_base& if (!Ch_front_var_stationnaire::initialiser(temps,inco)) return 0; - int dim=nb_comp(); - double x,y,z; // int nbsf=faces.nb_som_faces(); - int k; - DoubleTab& tab=valeurs(); + int dim=nb_comp(); + DoubleTab& tab_tab=valeurs(); + DoubleTabView tab = tab_tab.view_wo(); const Domaine_VF& domaine_vf=ref_cast(Domaine_VF,frontiere_dis().domaine_dis()); - const DoubleTab& xv=domaine_vf.xv(); + CDoubleTabView xv = domaine_vf.xv().view_ro(); const Front_VF& le_bord = ref_cast(Front_VF,frontiere_dis()); - int nb_faces_bord_tot=le_bord.nb_faces_tot(); - ToDo_Kokkos("critical"); - for (int ind_face=0; ind_face1) - y=xv(face,1); - if(dimension>2) - z=xv(face,2); - - for( k=0; k(nb_faces_bord_tot), KOKKOS_LAMBDA(const int i) + { + double x,y,z; + int threadId = fxyzk.acquire(); + + int face; + if (i < nb_faces) + face = premiere_face + i; + else + face = faces_virt[i - nb_faces]; + + x = y = z = 0; + x = xv(face, 0); + if(l_dimension > 1) + y = xv(face, 1); + if(l_dimension > 2) + z = xv(face, 2); + + fxyzk.setVar(0, x, threadId); + fxyzk.setVar(1, y, threadId); + fxyzk.setVar(2, z, threadId); + tab(i, k) = fxyzk.eval(threadId); + fxyzk.release(threadId); + //Cout << " x y z " << x << " " << y << " " << z << " " << tab(i,k) << finl; + }); + end_gpu_timer(__KERNEL_NAME__); } - tab.echange_espace_virtuel(); + tab_tab.echange_espace_virtuel(); return 1; } diff --git a/src/VDF/Champs/Champ_Face_VDF_implementation.cpp b/src/VDF/Champs/Champ_Face_VDF_implementation.cpp index e62dd7f992..513d4c6a6a 100644 --- a/src/VDF/Champs/Champ_Face_VDF_implementation.cpp +++ b/src/VDF/Champs/Champ_Face_VDF_implementation.cpp @@ -182,31 +182,37 @@ DoubleTab& Champ_Face_VDF_implementation::valeur_aux_faces_post_impl(const Domai return result; } -DoubleVect& Champ_Face_VDF_implementation::valeur_aux_elems_compo(const DoubleTab& positions, const IntVect& les_polys, DoubleVect& val, int ncomp) const +DoubleVect& Champ_Face_VDF_implementation::valeur_aux_elems_compo(const DoubleTab& tab_positions, const IntVect& tab_les_polys, DoubleVect& tab_val, int ncomp) const { - assert(val.size_totale() >= les_polys.size()); + assert(tab_val.size_totale() >= tab_les_polys.size()); const int D = Objet_U::dimension; - const DoubleTab& coord = domaine_vdf().domaine().coord_sommets(); - const IntTab& f_s = domaine_vdf().face_sommets(), &e_f = domaine_vdf().elem_faces(); - const DoubleTab& vals = le_champ().valeurs(); - int size = les_polys.size(); - ToDo_Kokkos("critical"); - for(int p = 0; p < size; p++) - { - int e = les_polys(p); - if (e<0) - { - val(p) = 0; - } - else - { - const double val1 = vals(e_f(e, ncomp)), val2 = vals(e_f(e, D + ncomp)); - const int som0 = f_s(e_f(e, ncomp), 0), som1 = f_s(e_f(e, D + ncomp), 0); - const double psi = (positions(p, ncomp) - coord(som0, ncomp)) / (coord(som1, ncomp) - coord(som0, ncomp)); - val(p) = interpolation(val1, val2, psi); - } - } - return val; + int size = tab_les_polys.size(); + + CDoubleTabView positions = tab_positions.view_ro(); + CIntArrView les_polys = tab_les_polys.view_ro(); + DoubleArrView val = tab_val.view_wo(); + + CDoubleTabView coord = domaine_vdf().domaine().coord_sommets().view_ro(); + CIntTabView f_s = domaine_vdf().face_sommets().view_ro(); + CIntTabView e_f = domaine_vdf().elem_faces().view_ro(); + CDoubleArrView vals = static_cast(le_champ().valeurs()).view_ro(); + Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),range_1D(0, size), KOKKOS_LAMBDA(const int p) + { + int e = les_polys(p); + if (e<0) + { + val(p) = 0; + } + else + { + const double val1 = vals(e_f(e, ncomp)), val2 = vals(e_f(e, D + ncomp)); + const int som0 = f_s(e_f(e, ncomp), 0), som1 = f_s(e_f(e, D + ncomp), 0); + const double psi = (positions(p, ncomp) - coord(som0, ncomp)) / (coord(som1, ncomp) - coord(som0, ncomp)); + val(p) = interpolation(val1, val2, psi); + } + }); + end_gpu_timer(__KERNEL_NAME__); + return tab_val; } double Champ_Face_VDF_implementation::valeur_a_elem_compo(const DoubleVect& position, int e, int d) const @@ -493,145 +499,148 @@ int Champ_Face_VDF_implementation::imprime_Face(Sortie& os, int ncomp) const return 1; } -DoubleTab& Champ_Face_VDF_implementation::trace(const Frontiere_dis_base& fr, const DoubleTab& y, DoubleTab& x,int distant) const +DoubleTab& Champ_Face_VDF_implementation::trace(const Frontiere_dis_base& fr, const DoubleTab& tab_y, DoubleTab& tab_x,int distant) const { assert(distant==0); const Front_VF& fr_vf=ref_cast(Front_VF, fr); const Domaine_VDF& zvdf=domaine_vdf(); - const IntVect& ori = zvdf.orientation(); - const IntTab& face_voisins = zvdf.face_voisins(); - const IntTab& elem_faces = zvdf.elem_faces(); - int elem1,elem2; - int face,i,f1,f2,f3,f4; + CIntArrView ori = zvdf.orientation().view_ro(); + CIntTabView face_voisins = zvdf.face_voisins().view_ro(); + CIntTabView elem_faces = zvdf.elem_faces().view_ro(); + int num_first_face = fr_vf.num_premiere_face(); int nb_faces = fr_vf.nb_faces(); - if (x.dimension(0)!=nb_faces) + if (tab_x.dimension(0)!=nb_faces) { Cerr << "The number of faces " << nb_faces << " on the remote boundary " << fr.le_nom() << finl; - Cerr << "does not match the number of faces " << x.dimension(0) << " on the local boundary." << finl; + Cerr << "does not match the number of faces " << tab_x.dimension(0) << " on the local boundary." << finl; Cerr << "Please, check if the boundary condition is not applied on wrong boundaries." << finl; Process::exit(); } // assert(x.dimension(1)==Objet_U::dimension); - if (x.dimension(1) == 1) + DoubleTabView x = tab_x.view_wo(); + CDoubleArrView y = static_cast(tab_y).view_ro(); + if (tab_x.dimension(1) == 1) { - ToDo_Kokkos("critical"); - for (i=0; idomaine(); - const IntTab& elem_faces = domaine_VDF.elem_faces(); - const DoubleVect& alpha_t = diffusivite_turbulente().valeurs(); + Domaine_VDF_View dom_view(domaine_VDF); + CIntTabView elem_faces = domaine_VDF.elem_faces().view_ro(); + const DoubleVect& tab_alpha_t = diffusivite_turbulente().valeurs(); + CDoubleArrView alpha_t = tab_alpha_t.view_ro(); bool is_concentration = (equation().que_suis_je().debute_par("Convection_Diffusion_Concentration") || equation().que_suis_je().debute_par("Convection_Diffusion_Espece")); - ArrOfInt numfa(2*dimension); - ToDo_Kokkos("critical"); - for (int elem = 0; elem < domaine_VDF.nb_elem(); elem++) + const int l_dimension = Objet_U::dimension; + ArrOfInt tab_numfa(2*l_dimension); + IntArrView numfa = tab_numfa.view_rw(); + int nb_elem = domaine_VDF.nb_elem(); + // Prefetch alpha_(elem) on host to avoid virtual function inside Kokkok lambda + DoubleArrView alpha_view; + for (int elem = 0; elem < nb_elem; elem++) { - // choix du facteur - double rcp = 1.; - if (!is_concentration) - { - const int Ccp = sub_type(Champ_Uniforme, mon_equation->milieu().capacite_calorifique()); - const int Cr = sub_type(Champ_Uniforme, mon_equation->milieu().masse_volumique()); - const DoubleTab& tab_Cp = mon_equation->milieu().capacite_calorifique().valeurs(); - const DoubleTab& tab_r = mon_equation->milieu().masse_volumique().valeurs(); - rcp = tab_r(Cr ? 0 : elem, 0) * tab_Cp(Ccp ? 0 : elem, 0); - } - - double moy = 0.; - for (int i = 0; i < 2 * dimension; i++) numfa[i] = elem_faces(elem, i); - - // XXX : E Saikali j'ai corrige pour multi inco parce que c'etait 1/dx et pas 1/dx^2 ... donc attention si ecart ! - // c'etait comme ca : for (int d = 0; d < dimension; d++) moy += 1. / (domaine_VDF.dist_face(numfa[d], numfa[dimension + d], d)); - for (int d = 0; d < dimension; d++) - { - const double hd = domaine_VDF.dist_face(numfa[d], numfa[dimension + d], d); - moy += 1. / (hd * hd); - } - const double alpha_local = (alpha_(elem) + alpha_t(elem)) / rcp * moy; - coef = std::max(coef, alpha_local); + alpha_view(elem) = alpha_(elem); + } + if (!is_concentration) + { + const int Ccp = sub_type(Champ_Uniforme, mon_equation->milieu().capacite_calorifique()); + const int Cr = sub_type(Champ_Uniforme, mon_equation->milieu().masse_volumique()); + CDoubleTabView tab_Cp = mon_equation->milieu().capacite_calorifique().valeurs().view_ro(); + CDoubleTabView tab_r = mon_equation->milieu().masse_volumique().valeurs().view_ro(); + Kokkos::parallel_reduce(start_gpu_timer(__KERNEL_NAME__), range_1D(0, nb_elem), KOKKOS_LAMBDA (const int elem, double& local_max) + { + double rcp = tab_r(Cr ? 0 : elem, 0) * tab_Cp(Ccp ? 0 : elem, 0); + + double moy = 0.; + for (int i = 0; i < 2 * dimension; i++) numfa[i] = elem_faces(elem, i); + + // XXX : E Saikali j'ai corrige pour multi inco parce que c'etait 1/dx et pas 1/dx^2 ... donc attention si ecart ! + // c'etait comme ca : for (int d = 0; d < dimension; d++) moy += 1. / (domaine_VDF.dist_face(numfa[d], numfa[dimension + d], d)); + for (int d = 0; d < dimension; d++) + { + const double hd = dom_view.dist_face(numfa[d], numfa[dimension + d], d); + moy += 1. / (hd * hd); + } + const double alpha_local = (alpha_view(elem) + alpha_t(elem)) / rcp * moy; + local_max = std::max(coef, alpha_local); + }, + Kokkos::Max(coef) + ); + end_gpu_timer(__KERNEL_NAME__); + } + else + { + Kokkos::parallel_reduce(start_gpu_timer(__KERNEL_NAME__), range_1D(0, nb_elem), KOKKOS_LAMBDA (const int elem, double& local_max) + { + double rcp = 1.0f; + + double moy = 0.; + for (int i = 0; i < 2 * dimension; i++) numfa[i] = elem_faces(elem, i); + + // XXX : E Saikali j'ai corrige pour multi inco parce que c'etait 1/dx et pas 1/dx^2 ... donc attention si ecart ! + // c'etait comme ca : for (int d = 0; d < dimension; d++) moy += 1. / (domaine_VDF.dist_face(numfa[d], numfa[dimension + d], d)); + for (int d = 0; d < dimension; d++) + { + const double hd = dom_view.dist_face(numfa[d], numfa[dimension + d], d); + moy += 1. / (hd * hd); + } + const double alpha_local = (alpha_view(elem) + alpha_t(elem)) / rcp * moy; + local_max = std::max(coef, alpha_local); + }, + Kokkos::Max(coef) + ); + end_gpu_timer(__KERNEL_NAME__); } coef = Process::mp_max(coef);