Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

清理backflow稳定化中对于dot sol的依赖关系 #155

Merged
merged 11 commits into from
Jul 29, 2024
Merged

清理backflow稳定化中对于dot sol的依赖关系 #155

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
914825a
remove dot_sol in NS local assembly
ju-liu Jul 19, 2024
7fa7ec9
Merge branch 'master' into sun-dev
ju-liu Jul 19, 2024
6ca5cf5
global assembly in NS remove the dot sol dependency
ju-liu Jul 19, 2024
1221150
modify the interface class accordingly
ju-liu Jul 19, 2024
c9f6ada
modify the call of BackFlow assembly routine
ju-liu Jul 19, 2024
d929388
ruc project clean the local assem routien
ju-liu Jul 19, 2024
d6494ce
glboal assem in ruc calling sequence adjustmment
ju-liu Jul 19, 2024
464071f
adjust the global assem body
ju-liu Jul 19, 2024
5a9d8f9
make use of vector3 to simplify code
ju-liu Jul 20, 2024
07a8cbc
use vector3 to simplify code
ju-liu Jul 20, 2024
6e5d05a
vector3 class simplifies the code
ju-liu Jul 20, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 1 addition & 3 deletions examples/ns/include/PGAssem_NS_FEM.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -168,16 +168,14 @@ class PGAssem_NS_FEM : public IPGAssem
const ALocal_EBC * const &ebc_part );

// Backflow integral on outlet surfaces
void BackFlow_G( const PDNSolution * const &dot_sol,
const PDNSolution * const &sol,
void BackFlow_G( const PDNSolution * const &sol,
IPLocAssem * const &lassem_ptr,
FEAElement * const &element_s,
const IQuadPts * const &quad_s,
const ALocal_NBC * const &nbc_part,
const ALocal_EBC * const &ebc_part );

void BackFlow_KG( const double &dt,
const PDNSolution * const &dot_sol,
const PDNSolution * const &sol,
IPLocAssem * const &lassem_ptr,
FEAElement * const &element_s,
Expand Down
2 changes: 0 additions & 2 deletions examples/ns/include/PLocAssem_VMS_NS_GenAlpha.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -122,7 +122,6 @@ class PLocAssem_VMS_NS_GenAlpha : public IPLocAssem
const IQuadPts * const &quad );

virtual void Assem_Residual_BackFlowStab(
const double * const &dot_sol,
const double * const &sol,
FEAElement * const &element,
const double * const &eleCtrlPts_x,
Expand All @@ -132,7 +131,6 @@ class PLocAssem_VMS_NS_GenAlpha : public IPLocAssem

virtual void Assem_Tangent_Residual_BackFlowStab(
const double &dt,
const double * const &dot_sol,
const double * const &sol,
FEAElement * const &element,
const double * const &eleCtrlPts_x,
Expand Down
50 changes: 18 additions & 32 deletions examples/ns/src/PGAssem_NS_FEM.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -319,7 +319,7 @@ void PGAssem_NS_FEM::Assem_residual(
delete [] row_index; row_index = nullptr;

// Backflow stabilization residual contribution
BackFlow_G( sol_a, sol_b, lassem_ptr, elements, quad_s, nbc_part, ebc_part );
BackFlow_G( sol_b, lassem_ptr, elements, quad_s, nbc_part, ebc_part );

// Resistance type boundary condition
NatBC_Resis_G( curr_time, dt, dot_sol_np1, sol_np1, lassem_ptr, elements, quad_s,
Expand Down Expand Up @@ -410,7 +410,7 @@ void PGAssem_NS_FEM::Assem_tangent_residual(
delete [] row_index; row_index = nullptr;

// Backflow stabilization residual & tangent contribution
BackFlow_KG( dt, sol_a, sol_b, lassem_ptr, elements, quad_s, nbc_part, ebc_part );
BackFlow_KG( dt, sol_b, lassem_ptr, elements, quad_s, nbc_part, ebc_part );

// Resistance type boundary condition
NatBC_Resis_KG( curr_time, dt, dot_sol_np1, sol_np1, lassem_ptr, elements, quad_s,
Expand Down Expand Up @@ -476,26 +476,22 @@ void PGAssem_NS_FEM::NatBC_G( const double &curr_time, const double &dt,
}

void PGAssem_NS_FEM::BackFlow_G(
const PDNSolution * const &dot_sol,
const PDNSolution * const &sol,
IPLocAssem * const &lassem_ptr,
FEAElement * const &element_s,
const IQuadPts * const &quad_s,
const ALocal_NBC * const &nbc_part,
const ALocal_EBC * const &ebc_part )
{
double * array_a = new double [nlgn * dof_sol];
double * array_b = new double [nlgn * dof_sol];
double * local_as = new double [dof_sol * snLocBas];
double * local_bs = new double [dof_sol * snLocBas];
double * array = new double [nlgn * dof_sol];
double * local = new double [snLocBas * dof_sol];
int * LSIEN = new int [snLocBas];
double * sctrl_x = new double [snLocBas];
double * sctrl_y = new double [snLocBas];
double * sctrl_z = new double [snLocBas];
PetscInt * srow_index = new PetscInt [dof_mat * snLocBas];

dot_sol->GetLocalArray( array_a );
sol->GetLocalArray( array_b );
sol->GetLocalArray( array );

for(int ebc_id = 0; ebc_id < num_ebc; ++ebc_id)
{
Expand All @@ -507,11 +503,10 @@ void PGAssem_NS_FEM::BackFlow_G(

ebc_part -> get_ctrlPts_xyz(ebc_id, ee, sctrl_x, sctrl_y, sctrl_z);

GetLocal(array_a, LSIEN, snLocBas, local_as);
GetLocal(array_b, LSIEN, snLocBas, local_bs);
GetLocal(array, LSIEN, snLocBas, local);

lassem_ptr->Assem_Residual_BackFlowStab( local_as, local_bs,
element_s, sctrl_x, sctrl_y, sctrl_z, quad_s);
lassem_ptr->Assem_Residual_BackFlowStab( local, element_s,
sctrl_x, sctrl_y, sctrl_z, quad_s);

for(int ii=0; ii<snLocBas; ++ii)
{
Expand All @@ -523,10 +518,8 @@ void PGAssem_NS_FEM::BackFlow_G(
}
}

delete [] array_a; array_a = nullptr;
delete [] array_b; array_b = nullptr;
delete [] local_as; local_as = nullptr;
delete [] local_bs; local_bs = nullptr;
delete [] array; array = nullptr;
delete [] local; local = nullptr;
delete [] LSIEN; LSIEN = nullptr;
delete [] sctrl_x; sctrl_x = nullptr;
delete [] sctrl_y; sctrl_y = nullptr;
Expand All @@ -535,26 +528,22 @@ void PGAssem_NS_FEM::BackFlow_G(
}

void PGAssem_NS_FEM::BackFlow_KG( const double &dt,
const PDNSolution * const &dot_sol,
const PDNSolution * const &sol,
IPLocAssem * const &lassem_ptr,
FEAElement * const &element_s,
const IQuadPts * const &quad_s,
const ALocal_NBC * const &nbc_part,
const ALocal_EBC * const &ebc_part )
{
double * array_a = new double [nlgn * dof_sol];
double * array_b = new double [nlgn * dof_sol];
double * local_as = new double [dof_sol * snLocBas];
double * local_bs = new double [dof_sol * snLocBas];
double * array = new double [nlgn * dof_sol];
double * local = new double [snLocBas * dof_sol];
int * LSIEN = new int [snLocBas];
double * sctrl_x = new double [snLocBas];
double * sctrl_y = new double [snLocBas];
double * sctrl_z = new double [snLocBas];
PetscInt * srow_index = new PetscInt [dof_mat * snLocBas];

dot_sol->GetLocalArray( array_a );
sol->GetLocalArray( array_b );
sol->GetLocalArray( array );

for(int ebc_id = 0; ebc_id < num_ebc; ++ebc_id)
{
Expand All @@ -566,11 +555,10 @@ void PGAssem_NS_FEM::BackFlow_KG( const double &dt,

ebc_part -> get_ctrlPts_xyz(ebc_id, ee, sctrl_x, sctrl_y, sctrl_z);

GetLocal(array_a, LSIEN, snLocBas, local_as);
GetLocal(array_b, LSIEN, snLocBas, local_bs);
GetLocal(array, LSIEN, snLocBas, local);

lassem_ptr->Assem_Tangent_Residual_BackFlowStab( dt, local_as, local_bs,
element_s, sctrl_x, sctrl_y, sctrl_z, quad_s);
lassem_ptr->Assem_Tangent_Residual_BackFlowStab( dt, local, element_s,
sctrl_x, sctrl_y, sctrl_z, quad_s);

for(int ii=0; ii<snLocBas; ++ii)
{
Expand All @@ -585,10 +573,8 @@ void PGAssem_NS_FEM::BackFlow_KG( const double &dt,
}
}

delete [] array_a; array_a = nullptr;
delete [] array_b; array_b = nullptr;
delete [] local_as; local_as = nullptr;
delete [] local_bs; local_bs = nullptr;
delete [] array; array = nullptr;
delete [] local; local = nullptr;
delete [] LSIEN; LSIEN = nullptr;
delete [] sctrl_x; sctrl_x = nullptr;
delete [] sctrl_y; sctrl_y = nullptr;
Expand Down
54 changes: 25 additions & 29 deletions examples/ns/src/PLocAssem_VMS_NS_GenAlpha.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -827,7 +827,6 @@ void PLocAssem_VMS_NS_GenAlpha::Assem_Residual_EBC_Resistance(
}

void PLocAssem_VMS_NS_GenAlpha::Assem_Residual_BackFlowStab(
const double * const &dot_sol,
const double * const &sol,
FEAElement * const &element,
const double * const &eleCtrlPts_x,
Expand All @@ -845,36 +844,34 @@ void PLocAssem_VMS_NS_GenAlpha::Assem_Residual_BackFlowStab(
{
const std::vector<double> R = element->get_R(qua);

double surface_area, factor;
double surface_area;

const Vector_3 n_out = element->get_2d_normal_out(qua, surface_area);

double u = 0.0, v = 0.0, w = 0.0;
Vector_3 velo(0.0, 0.0, 0.0);
for(int ii=0; ii<snLocBas; ++ii)
{
const int ii4 = ii * 4;
u += sol[ii4+1] * R[ii];
v += sol[ii4+2] * R[ii];
w += sol[ii4+3] * R[ii];
velo.x() += sol[ii4+1] * R[ii];
velo.y() += sol[ii4+2] * R[ii];
velo.z() += sol[ii4+3] * R[ii];
}

const double temp = u * n_out.x() + v * n_out.y() + w * n_out.z();
const double temp = Vec3::dot_product( velo, n_out );

if(temp < 0.0) factor = temp * rho0 * beta;
else factor = 0.0;
const double factor = temp < 0.0 ? temp * rho0 * beta : 0.0;

for(int A=0; A<snLocBas; ++A)
{
sur_Residual[4*A+1] -= surface_area * quad -> get_qw(qua) * R[A] * factor * u;
sur_Residual[4*A+2] -= surface_area * quad -> get_qw(qua) * R[A] * factor * v;
sur_Residual[4*A+3] -= surface_area * quad -> get_qw(qua) * R[A] * factor * w;
sur_Residual[4*A+1] -= surface_area * quad -> get_qw(qua) * R[A] * factor * velo.x();
sur_Residual[4*A+2] -= surface_area * quad -> get_qw(qua) * R[A] * factor * velo.y();
sur_Residual[4*A+3] -= surface_area * quad -> get_qw(qua) * R[A] * factor * velo.z();
}
}
}

void PLocAssem_VMS_NS_GenAlpha::Assem_Tangent_Residual_BackFlowStab(
const double &dt,
const double * const &dot_sol,
const double * const &sol,
FEAElement * const &element,
const double * const &eleCtrlPts_x,
Expand All @@ -894,32 +891,31 @@ void PLocAssem_VMS_NS_GenAlpha::Assem_Tangent_Residual_BackFlowStab(
{
const std::vector<double> R = element->get_R(qua);

double surface_area, factor;
double surface_area;

const Vector_3 n_out = element->get_2d_normal_out(qua, surface_area);

double u = 0.0, v = 0.0, w = 0.0;
Vector_3 velo(0.0, 0.0, 0.0);
for(int ii=0; ii<snLocBas; ++ii)
{
u += sol[ii*4+1] * R[ii];
v += sol[ii*4+2] * R[ii];
w += sol[ii*4+3] * R[ii];
velo.x() += sol[ii*4+1] * R[ii];
velo.y() += sol[ii*4+2] * R[ii];
velo.z() += sol[ii*4+3] * R[ii];
}

const double temp = u * n_out.x() + v * n_out.y() + w * n_out.z();
const double temp = Vec3::dot_product( velo, n_out );

if(temp < 0.0) factor = temp * rho0 * beta;
else factor = 0.0;
const double factor = temp < 0.0 ? temp * rho0 * beta : 0.0;

const double gwts = surface_area * quad -> get_qw(qua);

// snLocBas = 3 for linear tri element
// 6 for quadratic tri element
for(int A=0; A<snLocBas; ++A)
{
sur_Residual[4*A+1] -= gwts * R[A] * factor * u;
sur_Residual[4*A+2] -= gwts * R[A] * factor * v;
sur_Residual[4*A+3] -= gwts * R[A] * factor * w;
sur_Residual[4*A+1] -= gwts * R[A] * factor * velo.x();
sur_Residual[4*A+2] -= gwts * R[A] * factor * velo.y();
sur_Residual[4*A+3] -= gwts * R[A] * factor * velo.z();

for(int B=0; B<snLocBas; ++B)
{
Expand Down Expand Up @@ -953,15 +949,15 @@ double PLocAssem_VMS_NS_GenAlpha::get_flowrate( const double * const &sol,
double surface_area;
const Vector_3 n_out = element->get_2d_normal_out(qua, surface_area);

double u = 0.0, v = 0.0, w = 0.0;
Vector_3 velo(0.0, 0.0, 0.0);
for(int ii=0; ii<snLocBas; ++ii)
{
u += sol[ii*4+1] * R[ii];
v += sol[ii*4+2] * R[ii];
w += sol[ii*4+3] * R[ii];
velo.x() += sol[ii*4+1] * R[ii];
velo.y() += sol[ii*4+2] * R[ii];
velo.z() += sol[ii*4+3] * R[ii];
}

flrate += surface_area * quad->get_qw(qua) * ( u * n_out.x() + v * n_out.y() + w * n_out.z() );
flrate += surface_area * quad->get_qw(qua) * Vec3::dot_product( velo, n_out );
}

return flrate;
Expand Down
4 changes: 1 addition & 3 deletions examples/ruc_fsi/include/PGAssem_Tet_CMM_GenAlpha.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -213,8 +213,7 @@ class PGAssem_Tet_CMM_GenAlpha : public IPGAssem
const ALocal_EBC * const &ebc_part );

// Backflow integral on outlet surfaces
void BackFlow_G( const PDNSolution * const &dot_sol,
const PDNSolution * const &sol,
void BackFlow_G( const PDNSolution * const &sol,
IPLocAssem * const &lassem_ptr,
FEAElement * const &element_s,
const IQuadPts * const &quad_s,
Expand All @@ -223,7 +222,6 @@ class PGAssem_Tet_CMM_GenAlpha : public IPGAssem
const ALocal_EBC * const &ebc_part );

void BackFlow_KG( const double &dt,
const PDNSolution * const &dot_sol,
const PDNSolution * const &sol,
IPLocAssem * const &lassem_ptr,
FEAElement * const &element_s,
Expand Down
2 changes: 0 additions & 2 deletions examples/ruc_fsi/include/PLocAssem_Tet_CMM_GenAlpha.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -120,7 +120,6 @@ class PLocAssem_Tet_CMM_GenAlpha : public IPLocAssem
const IQuadPts * const &quad );

virtual void Assem_Residual_BackFlowStab(
const double * const &dot_sol,
const double * const &sol,
FEAElement * const &element,
const double * const &eleCtrlPts_x,
Expand All @@ -130,7 +129,6 @@ class PLocAssem_Tet_CMM_GenAlpha : public IPLocAssem

virtual void Assem_Tangent_Residual_BackFlowStab(
const double &dt,
const double * const &dot_sol,
const double * const &sol,
FEAElement * const &element,
const double * const &eleCtrlPts_x,
Expand Down
Loading