#include <godunovmethod.h>
Public Member Functions | |
GodunovMethod (OrthoMesh &mesh, Function3D &fInitU, const VecDouble &cPor, Function3D &fPrescribedU, FaceFluxFunction &flux, FixedValueCondition &fixedC, double CFL) | |
virtual void | iterateN (unsigned nSteps, double dt) |
Definition at line 10 of file godunovmethod.h.
GodunovMethod::GodunovMethod | ( | OrthoMesh & | mesh, | |
Function3D & | fInitU, | |||
const VecDouble & | cPor, | |||
Function3D & | fPrescribedU, | |||
FaceFluxFunction & | flux, | |||
FixedValueCondition & | fixedC, | |||
double | CFL | |||
) |
Definition at line 4 of file godunovmethod.cpp.
00005 :LaxFriedrichsForSystem(mesh,fInitU,cPor,fPrescribedU,flux,fixedC,CFL) 00006 { 00007 00008 }
void GodunovMethod::iterateN | ( | unsigned | nSteps, | |
double | dt | |||
) | [virtual] |
Reimplemented from LaxFriedrichsForSystem.
Definition at line 12 of file godunovmethod.cpp.
00013 { 00014 unsigned posCell,negCell; 00015 00016 //We iterate the method along faces, so 00017 //alloc the vectors to contain the solutions in neg 00018 //and pos cells of the faces. Remeber that the solution 00019 //in a cell is vector with a value for each component 00020 static VecDouble vBC(n_components),vFlux(n_components); 00021 static VecDoubleRef vQPos(n_components),vQNeg(n_components); 00022 00023 //For each time step 00024 for (unsigned count = 0; count < nSteps; count++) 00025 { 00026 //Get the face iterator and the number of faces. 00027 OrthoMesh::Cash_Face_It face = m_mesh.begin_cash_face(); 00028 OrthoMesh::Cash_Face_It endf = m_mesh.end_cash_face(); 00029 m_vbc.rewind(); 00030 bool hasbc; 00031 //m_solAnt receives the solution in time tn since 00032 //m_sol will contain the solution in tn+1. 00033 m_solAnt = m_sol; 00034 00035 //For each face 00036 for (;face!=endf;face++) 00037 { 00038 face->getAdjCellIndices(negCell,posCell); 00039 //Get the values of the previous solution in the right and left of the face 00040 hasbc=this->getValuesOfTheFaceCells(m_mesh,m_vbc,m_solAnt,face->index(),negCell,posCell,vQNeg,vQPos); 00041 00042 //printf("face %d : Values %g %g %g | %g %g %g\n",face->index(),vQNeg(0),vQNeg(1),vQNeg(2),vQPos(0),vQPos(1),vQPos(2)); 00043 00044 //Get the indices of the cells that contain the face. 00045 //If the face is at boundary, it has just one cell. 00046 //In this case the method getAdjCells() return posCell == negCell 00047 if (posCell == OrthoMesh::INVALID_INDEX) 00048 posCell=negCell; 00049 else if (negCell == OrthoMesh::INVALID_INDEX) 00050 negCell=posCell; 00051 00052 00053 00054 double posPor = getPorosity()(posCell); 00055 double negPor = getPorosity()(negCell); 00056 double mPor = (posPor + negPor)/2.0; 00057 00058 00059 00060 m_flux.exactFluxAtFace(vFlux,face,vQNeg,vQPos); 00061 00062 00063 //Now update the cell values component by component 00064 for (unsigned cmp=0;cmp<n_components;cmp++) 00065 { 00066 00067 //See if the face has prescribed boundary condition BC 00068 double flux; 00069 if (!hasbc) 00070 { 00071 /* 00072 The face doesnt have boundary condition BC 00073 proceed to the usual computation of the flux 00074 to calculate the amount of the mass passed through 00075 the face and divide it by the cell volume 00076 */ 00077 00078 flux = dt*face->areaPerCellVol() * vFlux(cmp) - mPor*(vQPos(cmp)-vQNeg(cmp))/6.0; 00079 } 00080 else //the face has BC, so dont add the stability term 00081 { 00082 flux = dt*face->areaPerCellVol() * vFlux(cmp); 00083 } 00084 00085 00086 00087 //printf("face %d = %g,Pos: %d,Neg %d, SwPos = %g,SwNeg = %g \n",faceIndex,flux,face->getPosCellIndex(),face->getNegCellIndex(),SwPos,SwNeg); 00088 if (face->hasPosCell()) 00089 m_sol(posCell,cmp) += + flux/posPor; 00090 00091 if (face->hasNegCell()) 00092 m_sol(negCell,cmp) += - flux/negPor; 00093 00094 } //end for each component 00095 } //end for each face 00096 00097 00098 } //end for each time step 00099 00100 00101 }