From e32dd556def4c442dd0a5e59ab44d0b6d207e4d0 Mon Sep 17 00:00:00 2001 From: Johan Mabille Date: Thu, 24 Sep 2020 17:09:59 +0200 Subject: [PATCH 1/3] Replaced initialization and obvious computation with STL algorithms --- proteus/mprans/RANS2P2D.h | 50 ++++++++++++++++++++------------------- 1 file changed, 26 insertions(+), 24 deletions(-) diff --git a/proteus/mprans/RANS2P2D.h b/proteus/mprans/RANS2P2D.h index a490ad70df..5c6405fd79 100644 --- a/proteus/mprans/RANS2P2D.h +++ b/proteus/mprans/RANS2P2D.h @@ -24,6 +24,18 @@ const double DM3=1.0;//1-point-wise divergence, 0-point-wise rate of volume cha const double inertial_term=1.0; namespace proteus { + + namespace detail + { + template + inline std::array make_array(const T& t) + { + std::array res; + std::fill_n(res.begin(), N, t); + return res; + } + } + template using GeneralizedFunctions = equivalent_polynomials::GeneralizedFunctions_mix; @@ -1682,35 +1694,22 @@ namespace proteus for(int eN=0;eN(0.0); + auto elementResidual_p_check = detail::make_array(0.0); + auto elementResidual_mesh = detail::make_array(0.0); + auto elementResidual_u = detail::make_array(0.0); + auto elementResidual_v = detail::make_array(0.0); + auto pelementResidual_u = detail::make_array(0.0); + auto pelementResidual_v = detail::make_array(0.0); + auto velocityErrorElement = detail::make_array(0.0); + double eps_rho,eps_mu; bool element_active=false; elementIsActive[eN]=false; const double* elementResidual_w(NULL); double mesh_volume_conservation_element=0.0, mesh_volume_conservation_element_weak=0.0; - for (int i=0;i 0) { for (int I=0;I Date: Mon, 28 Sep 2020 14:11:51 +0200 Subject: [PATCH 2/3] Use pyarray to compute data offset --- proteus/mprans/RANS2P2D.h | 521 +++++++++++++++++++------------------- 1 file changed, 256 insertions(+), 265 deletions(-) diff --git a/proteus/mprans/RANS2P2D.h b/proteus/mprans/RANS2P2D.h index 5c6405fd79..88d819a752 100644 --- a/proteus/mprans/RANS2P2D.h +++ b/proteus/mprans/RANS2P2D.h @@ -13,6 +13,8 @@ #include "xtensor-python/pyarray.hpp" #include "mpi.h" +#include + namespace py = pybind11; #define ZEROVEC {0.,0.} @@ -1714,10 +1716,10 @@ namespace proteus { for (int I=0;I 1.0e-8) + if (fabs(p_trial_ref(k, vi) - p_trial[vi]) > 1.0e-8) { for (int vj=0; vj < nDOF_trial_element; vj++) - std::cout<<"Trial "< 1.0e-8) @@ -2052,10 +2052,10 @@ namespace proteus prob=true; } //velocity - if (fabs(vel_trial_ref.data()[k*nDOF_v_trial_element + vi] - vel_trial[vi]) > 1.0e-8) + if (fabs(vel_trial_ref(k, vi) - vel_trial[vi]) > 1.0e-8) { for (int vj=0; vj < nDOF_v_trial_element; vj++) - std::cout<<"Trial "< 1.0e-8) @@ -2115,7 +2115,7 @@ namespace proteus #else for (int j=0;j 0) { - int ball_index=get_distance_to_ball(nParticles, ball_center.data(), ball_radius.data(),x,y,z,distance_to_solids.data()[eN_k]); + int ball_index=get_distance_to_ball(nParticles, ball_center.data(), ball_radius.data(),x,y,z,distance_to_solids(eN, k)); get_normal_to_ith_ball(nParticles, ball_center.data(), ball_radius.data(),ball_index,x,y,z,ball_n[0],ball_n[1]); } else @@ -2160,30 +2160,24 @@ namespace proteus //distance_to_solids is given in Prestep } if (nParticles > 0) - phi_solid.data()[eN_k] = distance_to_solids.data()[eN_k]; - const double particle_eps = particle_epsFact*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]); + phi_solid(eN, k) = distance_to_solids(eN, k); + const double particle_eps = particle_epsFact*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter(eN)); // //calculate pde coefficients at quadrature points // - const double H_s = gf_s.H(particle_eps,phi_solid.data()[eN_k]); - const double D_s = gf_s.D(particle_eps,phi_solid.data()[eN_k]); + const double H_s = gf_s.H(particle_eps,phi_solid(eN, k)); + const double D_s = gf_s.D(particle_eps,phi_solid(eN, k)); if ( nParticles == 0 || H_s != 0.0 || D_s != 0.0) { element_active=true; elementIsActive[eN]=true; } //save velocity at quadrature points for other models to use - double p_e = q_u_0.data()[eN_k] - p, - u_e = q_u_1.data()[eN_k] - u, - v_e = q_u_2.data()[eN_k] - v, + double p_e = q_u_0(eN ,k) - p, + u_e = q_u_1(eN ,k) - u, + v_e = q_u_2(eN, k) - v, velocity_e=sqrt(u_e*u_e + v_e*v_e); - - /* q_u_0.data()[eN_k] = p; */ - /* q_u_1.data()[eN_k] = u; */ - /* q_u_2.data()[eN_k] = v; */ - /* q_u_3.data()[eN_k] = 0.0; */ - double rho,nu; if (gf.useExact) { @@ -2225,18 +2219,18 @@ namespace proteus sigma, rho, nu, - elementDiameter.data()[eN], + elementDiameter(eN), smagorinskyConstant, turbulenceClosureModel, g.data(), useVF, - vf.data()[eN_k], - phi.data()[eN_k], - &normal_phi.data()[eN_k_nSpace], - kappa_phi.data()[eN_k], + vf(eN, k), + phi(eN, k), + &normal_phi(eN, k, 0), + kappa_phi(eN, k), //VRANS porosity, - phi_solid.data()[eN_k],//distance to solid + phi_solid(eN, k),//distance to solid p_old, u_old, v_old, @@ -2255,8 +2249,8 @@ namespace proteus v, w, LAG_LES, - q_eddy_viscosity.data()[eN_k], - q_eddy_viscosity_last.data()[eN_k], + q_eddy_viscosity(eN, k), + q_eddy_viscosity_last(eN, k), mom_u_acc, dmom_u_acc_u, mom_v_acc, @@ -2309,20 +2303,20 @@ namespace proteus dmom_w_ham_u, dmom_w_ham_v, dmom_w_ham_w, - forcex.data()[eN_k], - forcey.data()[eN_k], - forcez.data()[eN_k]); - q_rho.data()[eN_k] = rho; + forcex(eN, k), + forcey(eN, k), + forcez(eN, k)); + q_rho(eN, k) = rho; //VRANS - mass_source = q_mass_source.data()[eN_k]; + mass_source = q_mass_source(eN, k); //todo: decide if these should be lagged or not? updateDarcyForchheimerTerms_Ergun(NONCONSERVATIVE_FORM, /* linearDragFactor, */ /* nonlinearDragFactor, */ /* porosity, */ /* meanGrainSize, */ - q_dragAlpha.data()[eN_k], - q_dragBeta.data()[eN_k], + q_dragAlpha(eN, k), + q_dragBeta(eN, k), eps_rho, eps_mu, rho_0, @@ -2330,19 +2324,19 @@ namespace proteus rho_1, nu_1, useVF, - vf.data()[eN_k], - phi.data()[eN_k], + vf(eN, k), + phi(eN, k), u, v, w, - q_velocity_sge.data()[eN_k_nSpace+0], - q_velocity_sge.data()[eN_k_nSpace+1], - q_velocity_sge.data()[eN_k_nSpace+1],//dummy entry for 2D - eps_porous.data()[elementFlags.data()[eN]], - phi_porous.data()[eN_k], - q_velocity_porous.data()[eN_k_nSpace+0], - q_velocity_porous.data()[eN_k_nSpace+1], - q_velocity_porous.data()[eN_k_nSpace+1],//dummy entry for 2D + q_velocity_sge(eN, k, 0), + q_velocity_sge(eN, k, 1), + q_velocity_sge(eN, k, 1),//dummy entry for 2D + eps_porous(elementFlags(eN)), + phi_porous(eN, k), + q_velocity_porous(eN, k, 0), + q_velocity_porous(eN, k, 1), + q_velocity_porous(eN, k, 1),//dummy entry for 2D mom_u_source, mom_v_source, mom_w_source, @@ -2362,14 +2356,14 @@ namespace proteus rho_1, nu_1, useVF, - vf.data()[eN_k], - phi.data()[eN_k], + vf(eN, k), + phi(eN, k), porosity, c_mu, //mwf hack - q_turb_var_0.data()[eN_k], - q_turb_var_1.data()[eN_k], - &q_turb_var_grad_0.data()[eN_k_nSpace], - q_eddy_viscosity.data()[eN_k], + q_turb_var_0(eN, k), + q_turb_var_1(eN, k), + &q_turb_var_grad_0(eN, k, 0), + q_eddy_viscosity(eN, k), mom_uu_diff_ten, mom_vv_diff_ten, mom_ww_diff_ten, @@ -2417,17 +2411,17 @@ namespace proteus // //calculate time derivative at quadrature points // - if (q_dV_last.data()[eN_k] <= -100) - q_dV_last.data()[eN_k] = dV; - q_dV.data()[eN_k] = dV; + if (q_dV_last(eN, k) <= -100) + q_dV_last(eN, k) = dV; + q_dV(eN, k) = dV; ck.bdf(alphaBDF, - q_mom_u_acc_beta_bdf.data()[eN_k]*q_dV_last.data()[eN_k]/dV, + q_mom_u_acc_beta_bdf(eN, k)*q_dV_last(eN, k)/dV, mom_u_acc, dmom_u_acc_u, mom_u_acc_t, dmom_u_acc_u_t); ck.bdf(alphaBDF, - q_mom_v_acc_beta_bdf.data()[eN_k]*q_dV_last.data()[eN_k]/dV, + q_mom_v_acc_beta_bdf(eN, k)*q_dV_last(eN, k)/dV, mom_v_acc, dmom_v_acc_v, mom_v_acc_t, @@ -2443,20 +2437,20 @@ namespace proteus //calculate strong residual pdeResidual_p = ck.Advection_strong(dmass_adv_u,grad_u) + ck.Advection_strong(dmass_adv_v,grad_v) + - DM2*MOVING_DOMAIN*ck.Reaction_strong(alphaBDF*(dV-q_dV_last.data()[eN_k])/dV - div_mesh_velocity) + + DM2*MOVING_DOMAIN*ck.Reaction_strong(alphaBDF*(dV-q_dV_last(eN, k))/dV - div_mesh_velocity) + ck.Reaction_strong(mass_source); if (NONCONSERVATIVE_FORM > 0.0) { dmom_adv_sge[0] = 0.0; dmom_adv_sge[1] = 0.0; - dmom_ham_grad_sge[0] = inertial_term*dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+0] - MOVING_DOMAIN*xt); - dmom_ham_grad_sge[1] = inertial_term*dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+1] - MOVING_DOMAIN*yt); + dmom_ham_grad_sge[0] = inertial_term*dmom_u_acc_u*(q_velocity_sge(eN, k, 0) - MOVING_DOMAIN*xt); + dmom_ham_grad_sge[1] = inertial_term*dmom_u_acc_u*(q_velocity_sge(eN, k, 1) - MOVING_DOMAIN*yt); } else { - dmom_adv_sge[0] = inertial_term*dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+0] - MOVING_DOMAIN*xt); - dmom_adv_sge[1] = inertial_term*dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+1] - MOVING_DOMAIN*yt); + dmom_adv_sge[0] = inertial_term*dmom_u_acc_u*(q_velocity_sge(eN, k, 0) - MOVING_DOMAIN*xt); + dmom_adv_sge[1] = inertial_term*dmom_u_acc_u*(q_velocity_sge(eN, k, 1) - MOVING_DOMAIN*yt); dmom_ham_grad_sge[0] = 0.0; dmom_ham_grad_sge[1] = 0.0; } @@ -2482,7 +2476,7 @@ namespace proteus //add contributions from mass and source terms double tmpR=dmom_u_acc_u_t + dmom_u_source[0]; calculateSubgridError_tau(hFactor, - elementDiameter.data()[eN], + elementDiameter(eN), tmpR,//dmom_u_acc_u_t, dmom_u_acc_u, mv_tau,//dmom_adv_sge, @@ -2490,7 +2484,7 @@ namespace proteus dmom_u_ham_grad_p[0], tau_v0, tau_p0, - q_cfl.data()[eN_k]); + q_cfl(eN, k)); calculateSubgridError_tau(Ct_sge,Cd_sge, G,G_dd_G,tr_G, @@ -2500,7 +2494,7 @@ namespace proteus dmom_u_ham_grad_p[0], tau_v1, tau_p1, - q_cfl.data()[eN_k]); + q_cfl(eN, k)); tau_v = useMetrics*tau_v1+(1.0-useMetrics)*tau_v0; tau_p = useMetrics*tau_p1+(1.0-useMetrics)*tau_p0; @@ -2516,14 +2510,14 @@ namespace proteus subgridError_v, subgridError_w); // velocity used in adjoint (VMS or RBLES, with or without lagging the grid scale velocity) - dmom_adv_star[0] = inertial_term*dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+0] - MOVING_DOMAIN*xt + useRBLES*subgridError_u); - dmom_adv_star[1] = inertial_term*dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+1] - MOVING_DOMAIN*yt + useRBLES*subgridError_v); + dmom_adv_star[0] = inertial_term*dmom_u_acc_u*(q_velocity_sge(eN, k, 0) - MOVING_DOMAIN*xt + useRBLES*subgridError_u); + dmom_adv_star[1] = inertial_term*dmom_u_acc_u*(q_velocity_sge(eN, k, 1) - MOVING_DOMAIN*yt + useRBLES*subgridError_v); - mom_u_adv[0] += inertial_term*dmom_u_acc_u*(useRBLES*subgridError_u*q_velocity_sge.data()[eN_k_nSpace+0]); - mom_u_adv[1] += inertial_term*dmom_u_acc_u*(useRBLES*subgridError_v*q_velocity_sge.data()[eN_k_nSpace+0]); + mom_u_adv[0] += inertial_term*dmom_u_acc_u*(useRBLES*subgridError_u*q_velocity_sge(eN, k, 0)); + mom_u_adv[1] += inertial_term*dmom_u_acc_u*(useRBLES*subgridError_v*q_velocity_sge(eN, k, 0)); - mom_v_adv[0] += inertial_term*dmom_u_acc_u*(useRBLES*subgridError_u*q_velocity_sge.data()[eN_k_nSpace+1]); - mom_v_adv[1] += inertial_term*dmom_u_acc_u*(useRBLES*subgridError_v*q_velocity_sge.data()[eN_k_nSpace+1]); + mom_v_adv[0] += inertial_term*dmom_u_acc_u*(useRBLES*subgridError_u*q_velocity_sge(eN, k, 1)); + mom_v_adv[1] += inertial_term*dmom_u_acc_u*(useRBLES*subgridError_v*q_velocity_sge(eN, k, 1)); // adjoint times the test functions for (int i=0;i 0) { //cek todo, this needs to be fixed for not exact @@ -2575,10 +2569,10 @@ namespace proteus { for (int I=0;I= 0.0) + if (phi_solid(eN, k) >= 0.0) { p_LI = fmax(p_LI, fabs(p_e)); u_LI = fmax(u_LI, fabs(u_e)); @@ -2760,18 +2754,18 @@ namespace proteus { register int i_nSpace=i*nSpace; elementResidual_mesh[i] += H_s*H_f*(ck.Reaction_weak(1.0,p_test_dV[i]) - - ck.Reaction_weak(1.0,p_test_dV[i]*q_dV_last.data()[eN_k]/dV) - + ck.Reaction_weak(1.0,p_test_dV[i]*q_dV_last(eN, k)/dV) - ck.Advection_weak(mesh_vel,&p_grad_test_dV[i_nSpace])); elementResidual_p[i] += H_s*H_f*(ck.Advection_weak(mass_adv,&p_grad_test_dV[i_nSpace]) + ck.Hamiltonian_weak(mass_ham, p_test_dV[i]) + DM*MOVING_DOMAIN*(ck.Reaction_weak(alphaBDF*1.0,p_test_dV[i]) - - ck.Reaction_weak(alphaBDF*1.0,p_test_dV[i]*q_dV_last.data()[eN_k]/dV) - + ck.Reaction_weak(alphaBDF*1.0,p_test_dV[i]*q_dV_last(eN, k)/dV) - ck.Advection_weak(mesh_vel,&p_grad_test_dV[i_nSpace])) + ck.Reaction_weak(mass_source,p_test_dV[i])); if (nDOF_test_element == nDOF_v_test_element) { elementResidual_p[i] += - H_s*H_f*(PRESSURE_PROJECTION_STABILIZATION * ck.pressureProjection_weak(mom_uu_diff_ten[1], p, p_element_avg, p_test_ref.data()[k*nDOF_test_element+i], dV) + + H_s*H_f*(PRESSURE_PROJECTION_STABILIZATION * ck.pressureProjection_weak(mom_uu_diff_ten[1], p, p_element_avg, p_test_ref(k, i), dV) + (1 - PRESSURE_PROJECTION_STABILIZATION) * ck.SubgridError(subgridError_u,Lstar_u_p[i]) + (1 - PRESSURE_PROJECTION_STABILIZATION) * ck.SubgridError(subgridError_v,Lstar_v_p[i])); } @@ -2796,7 +2790,7 @@ namespace proteus ck.Reaction_weak(mom_u_source+NONCONSERVATIVE_FORM*dmom_u_acc_u*u*div_mesh_velocity,vel_test_dV[i]) + ck.Hamiltonian_weak(mom_u_ham,vel_test_dV[i]) + MOMENTUM_SGE*VELOCITY_SGE*ck.SubgridError(subgridError_u,Lstar_u_u[i]) + - ck.NumericalDiffusion(q_numDiff_u_last.data()[eN_k],grad_u,&vel_grad_test_dV[i_nSpace])); + ck.NumericalDiffusion(q_numDiff_u_last(eN, k),grad_u,&vel_grad_test_dV[i_nSpace])); elementResidual_v[i] += H_s*H_f*(ck.Mass_weak(mom_v_acc_t,vel_test_dV[i]) + ck.Advection_weak(mom_v_adv,&vel_grad_test_dV[i_nSpace]) + ck.Diffusion_weak(sdInfo_v_u_rowptr.data(),sdInfo_v_u_colind.data(),mom_vu_diff_ten,grad_u,&vel_grad_test_dV[i_nSpace]) + @@ -2804,7 +2798,7 @@ namespace proteus ck.Reaction_weak(mom_v_source+NONCONSERVATIVE_FORM*dmom_v_acc_v*v*div_mesh_velocity,vel_test_dV[i]) + ck.Hamiltonian_weak(mom_v_ham,vel_test_dV[i]) + MOMENTUM_SGE*VELOCITY_SGE*ck.SubgridError(subgridError_v,Lstar_v_v[i]) + - ck.NumericalDiffusion(q_numDiff_v_last.data()[eN_k],grad_v,&vel_grad_test_dV[i_nSpace])); + ck.NumericalDiffusion(q_numDiff_v_last(eN, k),grad_v,&vel_grad_test_dV[i_nSpace])); elementResidual_u[i] += H_s*H_f*MOMENTUM_SGE*PRESSURE_SGE*ck.SubgridError(subgridError_p,Lstar_p_u[i]); elementResidual_v[i] += H_s*H_f*MOMENTUM_SGE*PRESSURE_SGE*ck.SubgridError(subgridError_p,Lstar_p_v[i]); if (nParticles > 0)//solid boundary terms @@ -2818,12 +2812,12 @@ namespace proteus } }//i //estimate the numerical viscosity combining shock capturing and VMS/SUPG - numerical_viscosity.data()[eN_k] = q_numDiff_u_last.data()[eN_k] + MOMENTUM_SGE*VELOCITY_SGE*tau_v*(dmom_adv_star[0]*dmom_adv_star[0]+ + numerical_viscosity(eN, k) = q_numDiff_u_last(eN, k) + MOMENTUM_SGE*VELOCITY_SGE*tau_v*(dmom_adv_star[0]*dmom_adv_star[0]+ dmom_adv_star[1]*dmom_adv_star[1]); if (!elementIsActive[eN]) { - assert(std::fabs(gf_s.H(particle_eps,phi_solid.data()[eN_k])) == 0.0); - assert(std::fabs(gf_s.D(particle_eps,phi_solid.data()[eN_k])) == 0.0); + assert(std::fabs(gf_s.H(particle_eps,phi_solid(eN, k))) == 0.0); + assert(std::fabs(gf_s.D(particle_eps,phi_solid(eN, k))) == 0.0); } }//k }//fluid_phase @@ -2832,9 +2826,9 @@ namespace proteus { //compute indices and declare local storage register int eN_k = eN*nQuadraturePoints_element+k; - q_numDiff_u.data()[eN_k] = numDiffMax; - q_numDiff_v.data()[eN_k] = numDiffMax; - q_numDiff_w.data()[eN_k] = numDiffMax; + q_numDiff_u(eN, k) = numDiffMax; + q_numDiff_v(eN, k) = numDiffMax; + q_numDiff_w(eN, k) = numDiffMax; } #endif // @@ -2843,17 +2837,17 @@ namespace proteus for(int i=0;i DWp_Dn_jump, DW_Dn_jump; - register double gamma_cutfem=ghost_penalty_constant,gamma_cutfem_p=ghost_penalty_constant,h_cutfem=elementBoundaryDiameter.data()[*it]; - int eN_nDOF_v_trial_element = elementBoundaryElementsArray.data()[(*it)*2+0]*nDOF_v_trial_element; + register double gamma_cutfem=ghost_penalty_constant,gamma_cutfem_p=ghost_penalty_constant,h_cutfem=elementBoundaryDiameter(*it); + int eN = elementBoundaryElementsArray.data()[(*it)*2+0]; + int eN_nDOF_v_trial_element = eN*nDOF_v_trial_element; //See Massing Schott Wall 2018 //cek todo modify for two-fluids: rho_0 != rho_1 double norm_v=0.0; for (int i=0;i 0) { - get_distance_to_ball(nParticles, ball_center.data(), ball_radius.data(),x_ext,y_ext,z_ext,ebqe_phi_s.data()[ebNE_kb]); + get_distance_to_ball(nParticles, ball_center.data(), ball_radius.data(),x_ext,y_ext,z_ext,ebqe_phi_s(ebNE, kb)); } //else ebqe_phi_s.data()[ebNE_kb] is computed in Prestep const double particle_eps = particle_epsFact*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]); @@ -3356,19 +3350,19 @@ namespace proteus sigma, rho, nu, - elementDiameter.data()[eN], + elementDiameter(eN), smagorinskyConstant, turbulenceClosureModel, g.data(), useVF, - ebqe_vf_ext.data()[ebNE_kb], - ebqe_phi_ext.data()[ebNE_kb], - &ebqe_normal_phi_ext.data()[ebNE_kb_nSpace], - ebqe_kappa_phi_ext.data()[ebNE_kb], + ebqe_vf_ext(ebNE, kb), + ebqe_phi_ext(ebNE, kb), + &ebqe_normal_phi_ext(ebNE, kb, 0), + ebqe_kappa_phi_ext(ebNE, kb), //VRANS porosity_ext, // - ebqe_phi_s.data()[ebNE_kb], + ebqe_phi_s(ebNE, kb), p_old, u_old, v_old, @@ -3386,8 +3380,8 @@ namespace proteus v_ext, w_ext, LAG_LES, - ebqe_eddy_viscosity.data()[ebNE_kb], - ebqe_eddy_viscosity_last.data()[ebNE_kb], + ebqe_eddy_viscosity(ebNE, kb), + ebqe_eddy_viscosity_last(ebNE, kb), mom_u_acc_ext, dmom_u_acc_u_ext, mom_v_acc_ext, @@ -3453,19 +3447,19 @@ namespace proteus sigma, rho, nu, - elementDiameter.data()[eN], + elementDiameter(eN), smagorinskyConstant, turbulenceClosureModel, g.data(), useVF, - bc_ebqe_vf_ext.data()[ebNE_kb], - bc_ebqe_phi_ext.data()[ebNE_kb], - &ebqe_normal_phi_ext.data()[ebNE_kb_nSpace], - ebqe_kappa_phi_ext.data()[ebNE_kb], + bc_ebqe_vf_ext(ebNE, kb), + bc_ebqe_phi_ext(ebNE, kb), + &ebqe_normal_phi_ext(ebNE, kb, 0), + ebqe_kappa_phi_ext(ebNE, kb), //VRANS porosity_ext, // - ebqe_phi_s.data()[ebNE_kb], + ebqe_phi_s(ebNE, kb), p_old, u_old, v_old, @@ -3484,7 +3478,7 @@ namespace proteus bc_w_ext, LAG_LES, bc_eddy_viscosity_ext, - ebqe_eddy_viscosity_last.data()[ebNE_kb], + ebqe_eddy_viscosity_last(ebNE, kb), bc_mom_u_acc_ext, bc_dmom_u_acc_u_ext, bc_mom_v_acc_ext, @@ -3555,14 +3549,14 @@ namespace proteus rho_1, nu_1, useVF, - ebqe_vf_ext.data()[ebNE_kb], - ebqe_phi_ext.data()[ebNE_kb], + ebqe_vf_ext(ebNE, kb), + ebqe_phi_ext(ebNE, kb), porosity_ext, c_mu, //mwf hack - ebqe_turb_var_0.data()[ebNE_kb], - ebqe_turb_var_1.data()[ebNE_kb], + ebqe_turb_var_0(ebNE, kb), + ebqe_turb_var_1(ebNE, kb), turb_var_grad_0_dummy, //not needed - ebqe_eddy_viscosity.data()[ebNE_kb], + ebqe_eddy_viscosity(ebNE, kb), mom_uu_diff_ten_ext, mom_vv_diff_ten_ext, mom_ww_diff_ten_ext, @@ -3585,12 +3579,12 @@ namespace proteus rho_1, nu_1, useVF, - bc_ebqe_vf_ext.data()[ebNE_kb], - bc_ebqe_phi_ext.data()[ebNE_kb], + bc_ebqe_vf_ext(ebNE, kb), + bc_ebqe_phi_ext(ebNE, kb), porosity_ext, c_mu, //mwf hack - ebqe_turb_var_0.data()[ebNE_kb], - ebqe_turb_var_1.data()[ebNE_kb], + ebqe_turb_var_0(ebNE, kb), + ebqe_turb_var_1(ebNE, kb), turb_var_grad_0_dummy, //not needed bc_eddy_viscosity_ext, bc_mom_uu_diff_ten_ext, @@ -3653,16 +3647,16 @@ namespace proteus //calculate the numerical fluxes // ck.calculateGScale(G,normal,h_penalty); - penalty = useMetrics*C_b/h_penalty + (1.0-useMetrics)*ebqe_penalty_ext.data()[ebNE_kb]; + penalty = useMetrics*C_b/h_penalty + (1.0-useMetrics)*ebqe_penalty_ext(ebNE, kb); exteriorNumericalAdvectiveFlux(NONCONSERVATIVE_FORM, - isDOFBoundary_p.data()[ebNE_kb], - isDOFBoundary_u.data()[ebNE_kb], - isDOFBoundary_v.data()[ebNE_kb], - isDOFBoundary_w.data()[ebNE_kb], - isAdvectiveFluxBoundary_p.data()[ebNE_kb], - isAdvectiveFluxBoundary_u.data()[ebNE_kb], - isAdvectiveFluxBoundary_v.data()[ebNE_kb], - isAdvectiveFluxBoundary_w.data()[ebNE_kb], + isDOFBoundary_p(ebNE, kb), + isDOFBoundary_u(ebNE, kb), + isDOFBoundary_v(ebNE, kb), + isDOFBoundary_w(ebNE, kb), + isAdvectiveFluxBoundary_p(ebNE, kb), + isAdvectiveFluxBoundary_u(ebNE, kb), + isAdvectiveFluxBoundary_v(ebNE, kb), + isAdvectiveFluxBoundary_w(ebNE, kb), dmom_u_ham_grad_p_ext[0],//=1/rho, bc_dmom_u_ham_grad_p_ext[0],//=1/bc_rho, normal, @@ -3673,10 +3667,10 @@ namespace proteus bc_mom_u_adv_ext, bc_mom_v_adv_ext, bc_mom_w_adv_ext, - ebqe_bc_flux_mass_ext.data()[ebNE_kb]+MOVING_DOMAIN*(xt_ext*normal[0]+yt_ext*normal[1]),//BC is relative mass flux - ebqe_bc_flux_mom_u_adv_ext.data()[ebNE_kb], - ebqe_bc_flux_mom_v_adv_ext.data()[ebNE_kb], - ebqe_bc_flux_mom_w_adv_ext.data()[ebNE_kb], + ebqe_bc_flux_mass_ext(ebNE, kb)+MOVING_DOMAIN*(xt_ext*normal[0]+yt_ext*normal[1]),//BC is relative mass flux + ebqe_bc_flux_mom_u_adv_ext(ebNE, kb), + ebqe_bc_flux_mom_v_adv_ext(ebNE, kb), + ebqe_bc_flux_mom_w_adv_ext(ebNE, kb), p_ext, u_ext, v_ext, @@ -3704,30 +3698,30 @@ namespace proteus flux_mom_u_adv_ext, flux_mom_v_adv_ext, flux_mom_w_adv_ext, - &ebqe_velocity.data()[ebNE_kb_nSpace]); + &ebqe_velocity(ebNE, kb, 0)); for (int I=0;I nDOF_pressure) - nDOF_pressure=p_l2g.data()[eN_i]; + auto tmp = p_l2g(eN, i); + if (tmp > nDOF_pressure) + nDOF_pressure=tmp; } } nDOF_pressure +=1; assert(p_dof.shape(0) == nDOF_pressure); //std::cout<<"nDOF_pressure "< p_LI) From b04f9c3634ec4d0e63a60ddf3dd2d96caa68a6cf Mon Sep 17 00:00:00 2001 From: Johan Mabille Date: Tue, 29 Sep 2020 13:40:00 +0200 Subject: [PATCH 3/3] Code cleanup --- proteus/mprans/RANS2P2D.h | 39 ++++++++++----------------------------- 1 file changed, 10 insertions(+), 29 deletions(-) diff --git a/proteus/mprans/RANS2P2D.h b/proteus/mprans/RANS2P2D.h index 88d819a752..8eccc7ef9a 100644 --- a/proteus/mprans/RANS2P2D.h +++ b/proteus/mprans/RANS2P2D.h @@ -1966,18 +1966,10 @@ namespace proteus ck.gradTrialFromRef(&p_grad_trial_ref(k, 0, 0),jacInv,p_grad_trial); ck_v.gradTrialFromRef(&vel_grad_trial_ref(k, 0, 0),jacInv,vel_grad_trial); - for (int i=0; i < nDOF_trial_element; i++) - { - p_trial[i] = p_trial_ref(k, i); - p_grad_trial_ib[i*nSpace + 0] = p_grad_trial[i*nSpace+0]; - p_grad_trial_ib[i*nSpace + 1] = p_grad_trial[i*nSpace+1]; - } - for (int i=0; i < nDOF_v_trial_element; i++) - { - vel_trial[i] = vel_trial_ref(k, i); - vel_grad_trial_ib[i*nSpace + 0] = vel_grad_trial[i*nSpace+0]; - vel_grad_trial_ib[i*nSpace + 1] = vel_grad_trial[i*nSpace+1]; - } + std::copy_n(p_trial_ref.data() + k*nDOF_trial_element, nDOF_trial_element, p_trial); + std::copy(std::begin(p_grad_trial), std::end(p_grad_trial), std::begin(p_grad_trial_ib)); + std::copy_n(vel_trial_ref.data() + k*nDOF_v_trial_element, nDOF_v_trial_element, vel_trial); + std::copy(std::begin(vel_grad_trial), std::end(vel_grad_trial), std::begin(vel_grad_trial_ib)); if (icase == 0) { @@ -2836,7 +2828,6 @@ namespace proteus // for(int i=0;i(0.0); + auto elementResidual_p = detail::make_array(0.0); + auto elementResidual_u = detail::make_array(0.0); + auto elementResidual_v = detail::make_array(0.0); + register double eps_rho,eps_mu; const double* elementResidual_w(NULL); - for (int i=0;i