00001 #include "biphasicdiff.h"
00002 #include "hdf5orthowriter.h"
00003
00004
00005 BiphasicDiff::BiphasicDiff(OrthoMesh &mesh, const VecDouble &vK, const VecDouble &vPor,Function1D &fMobPc, Function1D &DfPc,LinearSolver &solver)
00006 :m_mesh(mesh),m_solver(solver),m_vK(vK),m_vPor(vPor),m_fMobPc(fMobPc),m_DfPc(DfPc)
00007 {
00008 m_B.reinit(mesh.numCells());
00009 buildSparsityPattern(mesh,m_sparsePattern);
00010 m_A.reinit(m_sparsePattern);
00011 m_vCoeff.reinit(mesh.numCells());
00012 }
00013
00014
00015 void BiphasicDiff::iterate(double dt)
00016 {
00017 ArrayOfVecDouble &transSol = getTransport().getSolutionAtCells();
00018 assert(transSol.vecs_size() == 1);
00019 VecDoubleRef vS(&(transSol(0,0)),transSol.size());
00020
00021
00022
00023 Index nCells = m_mesh.numCells();
00024 double cellVol=m_mesh.begin_cell()->volume();
00025
00026
00027 m_A=0.0;
00028 m_B=0.0;
00029
00030 for (Index cellIndex=0;cellIndex < nCells; cellIndex++)
00031 {
00032
00033
00034 double Saq=vS(cellIndex);;
00035 double por=m_vPor(cellIndex);
00036
00037 m_A.add(cellIndex,cellIndex,por*cellVol/dt);
00038 m_B(cellIndex)+=por*Saq/dt*cellVol;
00039
00040
00041 m_vCoeff(cellIndex)=-m_vK(cellIndex)*m_DfPc(Saq)*m_fMobPc(Saq);
00042 printf("Flux= %g %g %g\n",m_vK(cellIndex),m_DfPc(Saq),m_fMobPc(Saq));
00043
00044 }
00045
00046
00047 MixedHybridBase::assemblySystem(m_A,m_B,m_mesh,m_bc,m_vCoeff);
00048
00049
00050 m_solver.solve(m_A,vS,m_B);
00051
00052
00053
00054 }
00055
00056
00057 void BiphasicDiff::printOutput()
00058 {
00059
00060 ArrayOfVecDouble &transSol = getTransport().getSolutionAtCells();
00061 VecDoubleRef vS(&(transSol(0,0)),transSol.size());
00062 VecDouble Vq(m_mesh.numFaces());
00063 VecDouble Vv(m_mesh.numVertices());
00064 MixedHybridBase::getNormalVelocitiesFromPressure(Vq,m_mesh,m_bc,m_vCoeff,vS);
00065 HDF5OrthoWriter &hdf5 = HDF5OrthoWriter::getHDF5OrthoWriter();
00066 Matrix M(m_mesh.numVertices(),3);
00067 m_mesh.projectVelocitiesInFacesToVertices(M,Vq);
00068 M.copyColumn(Vv,2);
00069 hdf5.writeScalarField(Vv,"DiffVelZ");
00070 }
00071
