proteus  1.8.1
C/C++/Fortran libraries
VMS.cpp
Go to the documentation of this file.
1 #include "MeshAdaptPUMI.h"
2 #include <PCU.h>
3 #include <petscksp.h>
4 
5 #include <apf.h>
6 #include <apfMesh.h>
7 #include <apfShape.h>
8 #include <apfDynamicMatrix.h>
9 #include <apfNumbering.h>
10 
11 #include <iostream>
12 #include <fstream>
13 
14 //Global variables used to make it easier to pass these variables from MeshAdaptPUMIDrvr
15 extern int approx_order; //shape function order
16 extern int int_order; //integration order
17 extern double nu_0,nu_1,rho_0,rho_1;
18 double dt_err;
19 
20 struct Inputs{
21  apf::Vector3 vel_vect;
22  apf::Matrix3x3 gij;
23  apf::Matrix3x3 grad_vel;
24  apf::Vector3 grad_pres;
25  double visc_val;
26  double density;
27  int nsd;
28  apf::MeshElement* element;
29  apf::Element* pres_elem;
30  //apf::Element* visc_elem;
31  apf::Element* velo_elem;
32  //apf::Element* vof_elem;
33  apf::Element* velo_elem_old;
34  apf::Matrix3x3 KJ;
35  double* g;
36 };
37 
38 double get_nu_err(struct Inputs info);
39 apf::Vector3 getResidual(apf::Vector3 qpt,struct Inputs &info);
40 
41 
42 inline void getProps(double*rho,double*nu,double deltaT)
43 //Function used to transfer MeshAdaptPUMIDrvr variables into global variables
44 {
45  rho_0 = rho[0];
46  nu_0 = nu[0];
47  rho_1 = rho[1];
48  nu_1 = nu[1];
49  dt_err = deltaT;
50  return;
51 }
52 
53 apf::Matrix3x3 getKJ(int nsd)
54 {
55 //Jacobian correction matrix for proper metric tensor
56  if(nsd==2)
57  {
58  apf::Matrix3x3 KJtemp(1.0,0.5,0.0,0.5,1.0,0.0,0.0,0.0,1.0); //Jacobian correction matrix for proper metric tensor
59  return KJtemp;
60  }
61  else
62  {
63  double c1 = pow(2,1./3.); //1.259921049894873e+00;
64  double c2 = c1/2.0; //6.299605249474365e-01;
65  apf::Matrix3x3 KJtemp(c1,c2,c2,c2,c1,c2,c2,c2,c1);
66  return KJtemp;
67  }
68 }
69 
70 void MeshAdaptPUMIDrvr::get_VMS_error(double &total_error_out)
71 {
72  if(PCU_Comm_Self()==0)
73  std::cout<<"The beginning of the VMS\n";
74  getProps(rho,nu,delta_T);
75  approx_order = approximation_order;
76  int_order = integration_order;
77  nsd = m->getDimension();
78 
79  //***** Get Solution Fields First *****//
80  //apf::Field* voff = m->findField("vof");
81  //assert(voff);
82  apf::Field* velf = m->findField("velocity");
83  assert(velf);
84  apf::Field* pref = m->findField("p");
85  assert(pref);
86  apf::Field* velf_old;
87  if(m->findField("velocity_old")!=NULL)
88  velf_old = m->findField("velocity_old");
89  else{
90  velf_old = velf;
91  if(PCU_Comm_Self()==0)
92  std::cout<<"WARNING: old velocity field not found. Will proceed as if unsteady term is 0. \n";
93  dt_err = 1.0;
94  }
95  assert(velf_old);
96  //***** *****//
97  if(PCU_Comm_Self()==0)
98  std::cout<<"Got the solution fields\n";
99  //***** Compute the viscosity field *****//
100 /*
101  apf::Field* visc = getViscosityField(voff);
102  if(PCU_Comm_Self()==0)
103  std::cout<<"Got viscosity fields \n";
104 */
105  freeField(vmsErrH1);
106 
107  apf::Field* vmsErr = apf::createField(m,"VMSL2",apf::SCALAR,apf::getVoronoiShape(nsd,1));
108  //vmsErrH1 = apf::createField(m,"VMSH1",apf::SCALAR,apf::getVoronoiShape(nsd,1));
109  vmsErrH1 = apf::createField(m,"VMSH1",apf::SCALAR,apf::getVoronoiShape(nsd,1));
110 
111  if(PCU_Comm_Self()==0)
112  std::cout<<"Created the error fields\n";
113  //Start computing element quantities
114  int numqpt; //number of quadrature points
115  int nshl; //number of local shape functions
116  int elem_type; //what type of topology
117  double weight; //value container for the weight at each qpt
118  double Jdet;
119  //apf::EntityShape* elem_shape;
120  apf::Vector3 qpt; //container for quadrature points
121  apf::MeshElement* element;
122  apf::Element* visc_elem, *pres_elem,*velo_elem,*vof_elem,*velo_elem_old;
123  apf::Matrix3x3 J; //actual Jacobian matrix
124  apf::Matrix3x3 invJ; //inverse of Jacobian
125  apf::Matrix3x3 KJ = getKJ(nsd);
126  apf::Matrix3x3 gij; //actual Jacobian matrix
127 
128 
129  //apf::DynamicMatrix invJ_copy;
130  //apf::NewArray <apf::DynamicVector> shdrv;
131  //apf::NewArray <apf::DynamicVector> shgval_copy;
132 
133  apf::MeshIterator* iter = m->begin(nsd); //loop over elements
134  apf::MeshEntity* ent;
135 
136  int count = 0 ; //for debugging purposes
137  //Loop over elements and compute the VMS error in the L_2 norm
138  //double VMSerrTotalL2 = 0.0;
139  double VMSerrTotalH1 = 0.0;
140  while( (ent = m->iterate(iter)) ){
141 
142  element = apf::createMeshElement(m,ent);
143  pres_elem = apf::createElement(pref,element);
144  velo_elem = apf::createElement(velf,element);
145  velo_elem_old = apf::createElement(velf_old,element);
146  //visc_elem = apf::createElement(visc,element);
147  //vof_elem = apf::createElement(voff,element);
148 
149  double strongResidualTauL2;
150  double tau_m;
151  double areaCheck=0.0;
152  numqpt=apf::countIntPoints(element,int_order);
153 
154 /*
155  //Start the quadrature loop
156  for(int k=0;k<numqpt;k++){
157  apf::getIntPoint(element,int_order,k,qpt); //get a quadrature point and store in qpt
158  apf::getJacobian(element,qpt,J); //evaluate the Jacobian at the quadrature point
159  weight = apf::getIntWeight(element,int_order,k);
160  J = apf::transpose(J); //Is PUMI still defined in this way?
161  if(nsd==2)
162  J[2][2] = 1.0; //this is necessary to avoid singular matrix
163  invJ = invert(J);
164  Jdet=fabs(apf::getJacobianDeterminant(J,nsd));
165  gij = apf::transpose(invJ)*(KJ*invJ);
166 
167  double pressure = apf::getScalar(pres_elem,qpt);
168  apf::Vector3 grad_pres;
169  apf::getGrad(pres_elem,qpt,grad_pres);
170  double visc_val = apf::getScalar(visc_elem,qpt);
171  apf::Vector3 vel_vect;
172  apf::getVector(velo_elem,qpt,vel_vect);
173  apf::Matrix3x3 grad_vel;
174  apf::getVectorGrad(velo_elem,qpt,grad_vel);
175  grad_vel = apf::transpose(grad_vel);
176 
177  double C1 = 2.0; //constants for stabilization term
178  double C2 = 36.0;
179 
180  double stabTerm1 = 0.0;
181  double stabTerm2 = 0.0;
182  for(int i=0;i<nsd;i++){
183  for(int j=0;j<nsd;j++){
184  stabTerm1 += vel_vect[i]*gij[i][j]*vel_vect[j];
185  stabTerm2 += gij[i][j]*gij[i][j];
186  }
187  }
188  stabTerm2 = C2*visc_val*visc_val*stabTerm2;
189  tau_m = 1/sqrt(stabTerm1 + stabTerm2);
190 
191  //compute residual
192  apf::Vector3 tempConv;
193  apf::Vector3 tempDiff;
194  tempConv.zero();
195  tempDiff.zero();
196  for(int i=0;i<nsd;i++){
197  for(int j=0;j<nsd;j++){
198  tempConv[i] = tempConv[i] + vel_vect[j]*grad_vel[i][j];
199  }
200  }
201  double density = getMPvalue(apf::getScalar(vof_elem,qpt),rho_0,rho_1);
202  apf::Vector3 tempResidual = (tempConv + grad_pres/density);
203  double tempVal = tempResidual.getLength();
204  strongResidualTauL2 = tau_m*tau_m*tempVal*tempVal*weight*Jdet;
205 
206  } //end qpt loop
207 
208  double VMSerrL2 = sqrt(strongResidualTauL2);
209  apf::setScalar(vmsErr,ent,0,VMSerrL2);
210 */
211 
212  //H1 error compute nu_err at centroid and compute residual
213  qpt[0] = 0.0; qpt[1] = 0.0; qpt[2] = 0.0; //ensure it is initialized to 0
214  for(int i=0; i<nsd; i++)
215  qpt[i] = 1./(nsd+1.0);
216 
217  //double density = getMPvalue(apf::getScalar(vof_elem,qpt),rho_0,rho_1);
218  double density = rho[localNumber(ent)];
219  Inputs info;
220  info.element = element;
221  info.pres_elem = pres_elem;
222  //info.visc_elem = visc_elem;
223  info.velo_elem = velo_elem;
224  //info.vof_elem = vof_elem;
225  info.velo_elem_old = velo_elem_old;
226  info.KJ = KJ;
227  info.nsd = nsd;
228  info.density = density;
229  info.g = &g[0];
230 
231  info.visc_val = nu[localNumber(ent)];
232  apf::Vector3 tempResidual = getResidual(qpt,info);
233  double tempVal = tempResidual.getLength();
234  apf::getJacobian(element,qpt,J);
235  if(nsd==2)
236  J[2][2] = 1.0; //this is necessary to avoid singular matrix
237  invJ = invert(J);
238  Jdet=fabs(apf::getJacobianDeterminant(J,nsd));
239  gij = apf::transpose(invJ)*(KJ*invJ);
240  info.gij = gij;
241 
242  double nu_err = get_nu_err(info);
243  double VMSerrH1 = nu_err*tempVal*sqrt(apf::measure(m,ent));
244  //std::cout<<std::scientific<<std::setprecision(15)<<"H1 error for element "<<count<<" nu_err "<<nu_err<<" error "<<VMSerrH1<<std::endl;
245  apf::setScalar(vmsErrH1,ent,0,VMSerrH1);
246 
247  //VMSerrTotalL2 = VMSerrTotalL2+VMSerrL2*VMSerrL2;
248  VMSerrTotalH1 = VMSerrTotalH1+VMSerrH1*VMSerrH1;
249 
250  //apf::destroyElement(visc_elem);
251  apf::destroyElement(pres_elem);
252  apf::destroyElement(velo_elem);
253  apf::destroyElement(velo_elem_old);
254  //apf::destroyElement(vof_elem);
255  apf::destroyMeshElement(element);
256  count++;
257  } //end loop over elements
258 
259  //PCU_Add_Doubles(&VMSerrTotalL2,1);
260  PCU_Add_Doubles(&VMSerrTotalH1,1);
261  if(PCU_Comm_Self()==0)
262  std::cout<<std::scientific<<std::setprecision(15)<<"Total Error H1 "<<sqrt(VMSerrTotalH1)<<std::endl;
263  total_error = sqrt(VMSerrTotalH1);
264  total_error_out = sqrt(VMSerrTotalH1);
265  apf::destroyField(vmsErr);
266  //apf::destroyField(visc);
267 } //end function
268 
269 double get_nu_err(struct Inputs info){
270  double stabTerm1 = 0.0;
271  double stabTerm2 = 0.0;
272  for(int i=0;i<info.nsd;i++){
273  for(int j=0;j<info.nsd;j++){
274  stabTerm1 += info.vel_vect[i]*info.gij[i][j]*info.vel_vect[j];
275  stabTerm2 += info.gij[i][j]*info.gij[i][j];
276  }
277  }
278  double C_nu = 3.0;
279  stabTerm2 = C_nu*info.visc_val*info.visc_val*sqrt(stabTerm2);
280  double nu_err = 1.0/sqrt(info.visc_val*sqrt(stabTerm1) + stabTerm2);
281  return nu_err;
282 }
283 
284 apf::Vector3 getResidual(apf::Vector3 qpt,struct Inputs &info){
285  apf::Vector3 grad_pres;
286  apf::getGrad(info.pres_elem,qpt,grad_pres);
287  //double visc_val = apf::getScalar(info.visc_elem,qpt);
288  apf::Vector3 vel_vect;
289  apf::getVector(info.velo_elem,qpt,vel_vect);
290  apf::Vector3 vel_vect_old;
291  apf::getVector(info.velo_elem_old,qpt,vel_vect_old);
292  apf::Matrix3x3 grad_vel;
293  apf::getVectorGrad(info.velo_elem,qpt,grad_vel);
294  grad_vel = apf::transpose(grad_vel);
295 
296  apf::Vector3 tempConv;
297  apf::Vector3 tempDiff;
298  tempConv.zero();
299  tempDiff.zero();
300  for(int i=0;i<info.nsd;i++){
301  for(int j=0;j<info.nsd;j++){
302  tempConv[i] = tempConv[i] + vel_vect[j]*grad_vel[i][j];
303  }
304  tempConv[i] = tempConv[i] - info.g[i]; //body force contribution
305  }
306  apf::Vector3 tempResidual = (tempConv + grad_pres/info.density);
307  //acceleration term
308  tempResidual = tempResidual + (vel_vect-vel_vect_old)/dt_err;
309  //std::cout<<"What is the acceleration contribution? "<<(vel_vect-vel_vect_old)/dt_err<<std::endl;
310 
311  info.vel_vect = vel_vect;
312  info.grad_vel = grad_vel;
313  //info.visc_val = visc_val;
314  info.grad_pres = grad_pres;
315 
316  return tempResidual;
317 }
Inputs::pres_elem
apf::Element * pres_elem
Definition: VMS.cpp:29
MeshAdaptPUMIDrvr::total_error
double total_error
Definition: MeshAdaptPUMI.h:141
MeshAdaptPUMI.h
nu_0
double nu_0
Definition: ErrorResidualMethod.cpp:22
Inputs::g
double * g
Definition: VMS.cpp:35
MeshAdaptPUMIDrvr::nsd
int nsd
Definition: MeshAdaptPUMI.h:96
Inputs::KJ
apf::Matrix3x3 KJ
Definition: VMS.cpp:34
get_nu_err
double get_nu_err(struct Inputs info)
Definition: VMS.cpp:269
Inputs::grad_pres
apf::Vector3 grad_pres
Definition: VMS.cpp:24
MeshAdaptPUMIDrvr::localNumber
int localNumber(apf::MeshEntity *e)
Definition: MeshConverter.cpp:105
Inputs::visc_val
double visc_val
Definition: VMS.cpp:25
approx_order
int approx_order
Definition: ErrorResidualMethod.cpp:20
Inputs::element
apf::MeshElement * element
Definition: VMS.cpp:28
Inputs::grad_vel
apf::Matrix3x3 grad_vel
Definition: VMS.cpp:23
int_order
int int_order
Definition: ErrorResidualMethod.cpp:21
Inputs::velo_elem_old
apf::Element * velo_elem_old
Definition: VMS.cpp:33
MeshAdaptPUMIDrvr::integration_order
int integration_order
Definition: MeshAdaptPUMI.h:138
Inputs::gij
apf::Matrix3x3 gij
Definition: VMS.cpp:22
dt_err
double dt_err
Definition: VMS.cpp:18
MeshAdaptPUMIDrvr::approximation_order
int approximation_order
Definition: MeshAdaptPUMI.h:137
getProps
void getProps(double *rho, double *nu, double deltaT)
Definition: VMS.cpp:42
rho_1
double rho_1
Definition: VMS.cpp:17
rho_0
double rho_0
Definition: VMS.cpp:17
Inputs
Definition: VMS.cpp:20
Inputs::vel_vect
apf::Vector3 vel_vect
Definition: VMS.cpp:21
MeshAdaptPUMIDrvr::get_VMS_error
void get_VMS_error(double &total_error_out)
Definition: VMS.cpp:70
Inputs::velo_elem
apf::Element * velo_elem
Definition: VMS.cpp:31
getKJ
apf::Matrix3x3 getKJ(int nsd)
Definition: VMS.cpp:53
Inputs::density
double density
Definition: VMS.cpp:26
Inputs::nsd
int nsd
Definition: VMS.cpp:27
nu_1
double nu_1
Definition: VMS.cpp:17
getResidual
apf::Vector3 getResidual(apf::Vector3 qpt, struct Inputs &info)
Definition: VMS.cpp:284
Generated on Fri Jul 1 2022 11:17:46 for proteus by doxygen 1.8.20