00001 #include "fquadratic.h"
00002 #include "probefunction1d.h"
00003 #include "unittests.h"
00004 #include "numericmethods.h"
00005 #include <fstream>
00006 #include "constvelocitymodule.h"
00007 #include "fcompoundvector.h"
00008 #include "fplane.h"
00009 #include "flashco2h2o.h"
00010 #include "flashessential.h"
00011 #include "flash.h"
00012 #include "flashco2brine.h"
00013 #include "units.h"
00014 #include "flashsimpleblackoil.h"
00015 #include "flash.h"
00016 #define TEST(EXP)\
00017 printf("Testing " #EXP ":");\
00018 if (EXP)\
00019 printf("...OK\n");\
00020 else\
00021 {printf("...ERROR\nQUITTING......\n");abort();}
00022
00023
00024
00025 void UnitTests::testProbeFunction()
00026 {
00027 double min,max;
00028 printf("Testing ProbeFunction1D module.....\n");
00029 ProbeFunction1D fp(*new FQuadratic(1,0,0),0,-1,1,10000);
00030 fp.print();
00031 fp.getMinMaxValues(-0.5,-0.2,min,max);
00032 double tol=1.e-08;
00033 assert(NumericMethods::a_equal(max,0.25,tol));
00034 assert(NumericMethods::a_equal(min,0.04,tol));
00035
00036 fp.getMinMaxValues(-0.1,0.2,min,max);
00037 assert(NumericMethods::a_equal(max,0.04,tol));
00038 assert(NumericMethods::a_equal(min,0,tol));
00039
00040 fp.getMinMaxValues(0.1,0.2,min,max);
00041 assert(NumericMethods::a_equal(max,0.04,tol));
00042 assert(NumericMethods::a_equal(min,0.01,tol));
00043 printf("done\n");
00044
00045 }
00046
00047
00048
00049
00050
00051
00052 void UnitTests::plotSystemFluxFunction(OrthoMesh::Face_It &face,unsigned faceIndex,FaceFluxFunction &flux,double c1I,double c1E,unsigned cmp1,double c2,unsigned nPoints,std::string strFileName)
00053 {
00054
00055 std::ofstream fout(strFileName.c_str());
00056 VecDouble vQ1(flux.numComponents());
00057 VecDouble vQ2(flux.numComponents());
00058 VecDouble vOut(flux.numComponents());
00059
00060 vQ1 = c2;
00061 vQ2 = c2;
00062
00063
00064 double dx = (c1E-c1I)/nPoints;
00065 double x=c1I;
00066 unsigned cell1,cell2;
00067 face->getAdjCellIndices(cell1,cell2);
00068 flux.fluxAtFace(vOut,face,vQ1,vQ2);
00069 for (unsigned i=0;i<nPoints;i++)
00070 {
00071 vQ1(cmp1) = x;
00072 vQ2(cmp1) = x;
00073 fout << x << " " << vOut(cmp1) << std::endl;
00074 x+=dx;
00075
00076 }
00077 fout.close();
00078 }
00079
00080
00081
00082
00083 void UnitTests::testSystemFluxFunction()
00084 {
00085 app.printBeginSection("Testing the flux function of the system",std::cout);
00086 OrthoMesh &mesh = app.getMesh();
00087 FaceFluxFunction *flux = app.getFluxFunction();
00088
00089 ConstVelocityModule mod(mesh,*(new FCompoundVector(new FPlane(0,0,0,1),new FPlane(0,0,0,1),new FPlane(0,0,0,1))));
00090 flux->updateDynamicData(mod);
00091
00092 OrthoMesh::Face_It face = mesh.begin_face();
00093 OrthoMesh::Face_It endf = mesh.end_face();
00094 for (;face != endf;face++)
00095 if (!face->at_boundary() && face->getNormalNonZeroComponent() == 0)
00096 {
00097 printf("Printing flux in X direction: file fluxX\n");
00098 plotSystemFluxFunction(face,face->index(),*flux,0,1,0,0,1000,"fluxXC1_0.tst");
00099 plotSystemFluxFunction(face,face->index(),*flux,0,1,0,0.2,1000,"fluxXC1_0.2.tst");
00100 plotSystemFluxFunction(face,face->index(),*flux,0,1,0,0.5,1000,"fluxXC1_0.5.tst");
00101 plotSystemFluxFunction(face,face->index(),*flux,0,1,0,0.7,1000,"fluxXC1_0.7.tst");
00102 plotSystemFluxFunction(face,face->index(),*flux,0,1,0,1,1000,"fluxXC1_1.tst");
00103
00104
00105 plotSystemFluxFunction(face,face->index(),*flux,0,1,1,0,1000,"fluxXC2_0.tst");
00106 plotSystemFluxFunction(face,face->index(),*flux,0,1,1,0.2,1000,"fluxXC2_0.2.tst");
00107 plotSystemFluxFunction(face,face->index(),*flux,0,1,1,0.5,1000,"fluxXC2_0.5.tst");
00108 plotSystemFluxFunction(face,face->index(),*flux,0,1,1,0.7,1000,"fluxXC2_0.7.tst");
00109 plotSystemFluxFunction(face,face->index(),*flux,0,1,1,1,1000,"fluxXC2_1.tst");
00110 break;
00111 }
00112
00113 face=mesh.begin_face();
00114 endf=mesh.end_face();
00115 for (;face != endf;face++)
00116 if (!face->at_boundary() && face->getNormalNonZeroComponent() == 2)
00117 {
00118 printf("Printing flux in Z direction:\nComponent 1:\n");
00119 plotSystemFluxFunction(face,face->index(),*flux,0,1,0,0,1000,"fluxZC1_0.tst");
00120 plotSystemFluxFunction(face,face->index(),*flux,0,1,0,0.2,1000,"fluxZC1_0.2.tst");
00121 plotSystemFluxFunction(face,face->index(),*flux,0,1,0,0.5,1000,"fluxZC1_0.5.tst");
00122 plotSystemFluxFunction(face,face->index(),*flux,0,1,0,0.7,1000,"fluxZC1_0.7.tst");
00123 plotSystemFluxFunction(face,face->index(),*flux,0,1,0,1,1000,"fluxZC1_1.tst");
00124
00125 printf("Printing flux in Z direction:\nComponent 2:\n");
00126 plotSystemFluxFunction(face,face->index(),*flux,0,1,1,0,1000,"fluxZC2_0.tst");
00127 plotSystemFluxFunction(face,face->index(),*flux,0,1,1,0.2,1000,"fluxZC2_0.2.tst");
00128 plotSystemFluxFunction(face,face->index(),*flux,0,1,1,0.5,1000,"fluxZC2_0.5.tst");
00129 plotSystemFluxFunction(face,face->index(),*flux,0,1,1,0.7,1000,"fluxZC2_0.7.tst");
00130 plotSystemFluxFunction(face,face->index(),*flux,0,1,1,1,1000,"fluxZC2_1.tst");
00131
00132 printf("Printing Max Char Velocity in Z direction:\nComponent 1:\n");
00133 plotSystemFluxMaxLocalChar(face,face->index(),*flux,0,1,0.5,0.5,0,1000,"fluxZdC1_0.5.tst");
00134 plotSystemFluxMaxLocalChar(face,face->index(),*flux,0,1,0.0,0.0,0,1000,"fluxZdC1_0.tst");
00135 plotSystemFluxMaxLocalChar(face,face->index(),*flux,0,0.5,0.0,0.0,1,10,"fluxZdC2_0.tst");
00136 plotSystemFluxMaxLocalChar(face,face->index(),*flux,0,0.5,0.5,0.5,1,10,"fluxZdC2_0.5.tst");
00137 plotSystemFluxMaxLocalChar(face,face->index(),*flux,0,0.5,0.9,0.9,1,10,"fluxZdC2_0.9.tst");
00138 plotSystemFluxMaxLocalChar(face,face->index(),*flux,0,0.5,0.9,0.9,1,10,"fluxZdC2_0.9.tst");
00139 plotSystemFluxMaxLocalChar(face,face->index(),*flux,0,0.5,0.1,0.1,1,10,"fluxZdC2_0-0.1.tst");
00140
00141
00142
00143 break;
00144
00145
00146
00147 }
00148 }
00149
00150
00151 void UnitTests::execute()
00152 {
00153 unsigned resp;
00154 app.printBeginSection("Unit Tests",std::cout);
00155 Point3D P0(0,0,0);
00156 Point3D P1(1,1,1);
00157 OrthoMesh &mesh= app.getMesh();
00158
00159 while (1)
00160 {
00161 printf("\n\n=========================\n");
00162 printf("Welcome to the unit tests\n\
00163 =========================\n\n\
00164 Please choose the module to test:\n\
00165 \n\
00166 1) Mesh\n\
00167 2) FlashData\n\
00168 3) FlashCO2Brine\n\
00169 4) FlashSimpleBlackOil\n\
00170 5) Allan Flash\n\
00171 6) Print Mobilities and Fractional Flux\n\
00172 7) Print Capillary Pressure function\n\
00173 8) Flash CO2 Brine test\n\
00174 0) Quit\n\
00175 \n\
00176 Resp:");
00177 std::cin >> resp;
00178 switch(resp)
00179 {
00180 case 0:
00181 return;
00182 case 1:
00183 testMesh(mesh);
00184 break;
00185 case 2:
00186 testFlashData();
00187 break;
00188 case 3:
00189 testFlashCO2Brine();
00190 break;
00191 case 4:
00192 testSimpleBlackOilFlash();
00193 break;
00194 case 5:
00195 testAllanFlash();
00196 break;
00197 case 6:
00198 plotMobilities();
00199 break;
00200 case 7:
00201 plotPc();
00202 break;
00203 case 8:
00204 calculateFlashCO2Brine();
00205 default:
00206
00207 printf("Oooops Invalid Option....Exiting!!\n");
00208 return;
00209 }
00210 }
00211
00212
00213 }
00214
00215
00216 void UnitTests::plotSystemFluxMaxLocalChar(OrthoMesh::Face_It &face,unsigned faceIndex,FaceFluxFunction &flux,double c1I,double c1E,double c2I,double c2E,unsigned cmp,unsigned nPoints,std::string strFileName)
00217 {
00218 std::ofstream fout(strFileName.c_str());
00219 VecDouble vQ1(flux.numComponents());
00220 VecDouble vQ2(flux.numComponents());
00221 VecDouble vOut(flux.numComponents());
00222 vQ1 = c2I;
00223 vQ2 = c2E;
00224
00225
00226 double dx = (c1E-c1I)/nPoints;
00227 double x=c1I;
00228 unsigned cell1,cell2;
00229 face->getAdjCellIndices(cell1,cell2);
00230 double a=flux.maxLocalCharVelocity(face,vQ1,vQ2);
00231 for (unsigned i=1;i<nPoints;i++)
00232 {
00233 vQ1(cmp) = x;
00234 vQ2(cmp) = x+dx;
00235 fout << x << " " << a << std::endl;
00236 x+=dx;
00237
00238 }
00239 fout.close();
00240
00241 }
00242
00243
00244 void UnitTests::testDelIndexLst()
00245 {
00246 printf("Testing DelIndexLst::adjust() method: ........");
00247 DelIndexLst lst;
00248 lst.addIndex(5);
00249 lst.addIndex(4);
00250 lst.addIndex(13);
00251 lst.addIndex(7);
00252 lst.close();
00253 DelIndexLst::Iterator it=lst.begin();
00254
00255 if (lst.adjust(2,it) != 2 ||
00256 lst.adjust(0,it) != 0 ||
00257 lst.adjust(8,it) != 5 ||
00258 lst.adjust(6,it) != 4 ||
00259 lst.adjust(22,it) != 18 ||
00260 lst.adjust(0,it) != 0 ||
00261 lst.adjust(22,it) != 18)
00262 printf("Error\n");
00263 else
00264 printf("Ok");
00265
00266
00267
00268 }
00269
00270
00271 void UnitTests::testFaceCellMap(OrthoMesh &mesh,OrthoMesh::Cell_It &cell,FaceDirection3D dir)
00272 {
00273 unsigned findex = cell->face_index(dir);
00274 OrthoMesh::Face_It face = mesh.get_face_from_raw_index(cell->face_raw_index(dir));
00275 OrthoMesh::Face_It face2 = mesh.get_face(cell->face_index(dir));
00276 if (findex != OrthoMesh::INVALID_INDEX)
00277 {
00278 if (face != cell->face(dir))
00279 printf("Error of consistence, cell->face() != get_face_from_raw_index(cell->face_raw_index()) at cell %u",cell->index());
00280 if (findex != face->index())
00281 printf("Error of consistence, cell->face_index() == %u and cell->face()->index() == %u\n",cell->index(),face->index());
00282 if (findex != face2->index())
00283 printf("Error of consistence, cell->face_index() == %u and OrthoMesh::get_face(cell->face_index()) == %u\n",cell->index(),face->index());
00284
00285 if (dir == LEFT_FACE || dir == FRONT_FACE || dir == BOTTOM_FACE)
00286 {
00287 if (cell->index() != mesh.get_face(findex)->getNegCell())
00288 {
00289 printf("Error checking consistency beetween cell %u and face %u\n",cell->index(),findex);
00290 }
00291 }
00292 else
00293 {
00294 if (cell->index() != mesh.get_face(findex)->getPosCell())
00295 {
00296 printf("Error checking consistency beetween cell %u and face %u\n",cell->index(),findex);
00297 }
00298 }
00299 }
00300 }
00301
00302 void testCellIndexInSequence(OrthoMesh &mesh)
00303 {
00304 OrthoMesh::Cell_It cell = mesh.begin_cell();
00305 OrthoMesh::Cell_It endc = mesh.end_cell();
00306 unsigned index=0;
00307 printf("Testing Cell Indices continuity...");
00308 for (;cell!=endc;cell++,index++)
00309 {
00310 if (cell->index() != index)
00311 {
00312 printf("Error:\nIndice is not continuous in cell %u\n",cell->index());
00313 }
00314 }
00315 printf("Ok.\n");
00316 }
00317
00318
00319 void testFaceIndexInSequence(OrthoMesh &mesh)
00320 {
00321 OrthoMesh::Face_It cell = mesh.begin_face();
00322 OrthoMesh::Face_It endc = mesh.end_face();
00323 unsigned index=0;
00324 printf("Testing Face Indices continuity...");
00325 for (;cell!=endc;cell++,index++)
00326 {
00327 if (cell->index() != index)
00328 {
00329 printf("Error:\nIndice is not continuous in face %u\n",cell->index());
00330 }
00331 }
00332 printf("Ok.\n");
00333 }
00334
00335
00340 void UnitTests::testMesh(OrthoMesh &mesh)
00341 {
00342 testCellIndexInSequence(mesh);
00343 testFaceIndexInSequence(mesh);
00344
00345 printf("Printing mesh info to file mesh.dat.........");
00346 mesh.print("mesh.dat");
00347 printf("Ok\n" );
00348
00349 OrthoMesh::Cell_It cell = mesh.begin_cell();
00350 OrthoMesh::Cell_It endc = mesh.end_cell();
00351 printf("Testing FACE X CELL Mapping consistency......");
00352 for(;cell!=endc;cell++)
00353 {
00354 testFaceCellMap(mesh,cell,LEFT_FACE);
00355 testFaceCellMap(mesh,cell,RIGHT_FACE);
00356 testFaceCellMap(mesh,cell,BOTTOM_FACE);
00357 testFaceCellMap(mesh,cell,UP_FACE);
00358 testFaceCellMap(mesh,cell,FRONT_FACE);
00359 testFaceCellMap(mesh,cell,BACK_FACE);
00360 }
00361 printf("Ok.\n");
00362
00363 printf("Testing Inner Face Advance Index....");
00364
00365 OrthoMesh::Raw_Face_It faceIt = mesh.begin_raw_face();
00366 OrthoMesh::Raw_Face_It endFace = mesh.end_raw_face();
00367 OrthoMesh::Inner_Raw_Face_It iIt = mesh.begin_inner_raw_face();
00368 OrthoMesh::Inner_Raw_Face_It iEnd = mesh.end_inner_raw_face();
00369
00370 for (;iIt != iEnd;iIt++)
00371 {
00372 faceIt++;
00373 while (faceIt->at_boundary())
00374 faceIt++;
00375
00376 if (iIt->getI() == faceIt->getI() &&
00377 iIt->getJ() == faceIt->getJ() &&
00378 iIt->getK() == faceIt->getK() &&
00379 iIt->index() == faceIt->index())
00380 {
00381
00382
00383
00384
00385 }
00386 else
00387 {
00388 printf("Found inconsistency between\n");
00389 printf("\tInner Face: %d,%d,%d normal %d index: %d\n",iIt->getI(),iIt->getJ(),iIt->getK(),iIt->getNormalOrientation(),iIt->index());
00390 printf("\t Face: %d,%d,%d normal %d index: %d\n",faceIt->getI(),faceIt->getJ(),faceIt->getK(),faceIt->getNormalOrientation(),faceIt->index());
00391 abort();
00392 }
00393 }
00394 printf("Ok\n");
00395
00396 printf("Testing the get_face(index) function...." );
00397 {
00398 unsigned c1,c2,c3,c4;
00399 OrthoMesh::Face_It faceIt = mesh.begin_face();
00400 OrthoMesh::Face_It endFace = mesh.end_face();
00401 for(;faceIt!=endFace;faceIt++)
00402 {
00403
00404 OrthoMesh::Face_It fc = mesh.get_face(faceIt->index());
00405 if (!(fc->getI() == faceIt->getI() &&
00406 fc->getJ() == faceIt->getJ() &&
00407 fc->getK() == faceIt->getK() &&
00408 fc->index() == faceIt->index()))
00409 {
00410 printf("Test Failed begin_face(%u) != face\n",faceIt->index());
00411 abort();
00412 }
00413 fc->getAdjCellIndices(c1,c2);
00414 faceIt->getAdjCellIndices(c3,c4);
00415 if (c1 != c3 || c2 != c4)
00416 {
00417 printf("Test Failed begin_face(%u) != face\n",faceIt->index());
00418 abort();
00419 }
00420 }
00421 }
00422 printf("Ok\n");
00423
00424 }
00425
00426
00427 void UnitTests::testAllanFlash()
00428 {
00429 double P = 275.0;
00430 double P_MPa = P/10.0;
00431 double T = 353.0;
00432 double R = 8.314472;
00433 double KD = 497.0;
00434 double vn = 32.1;
00435
00436 double xCO2Henry = P_MPa/KD * exp(-vn * P_MPa/(R * T));
00437
00438 for(double zCO2 = 0.0; zCO2 <= 1.0; zCO2 += 0.01)
00439 {
00440 Flash::BrineCO2Data data;
00441
00442 Flash::ProgressiveFlashBrineCO2(data, T, P, 0.0, 1.0 - zCO2, zCO2);
00443
00444 printf("%g %g\n", zCO2, data.xCO2);
00445 }
00446
00447 printf("\n%g\n", xCO2Henry);
00448 }
00449
00450
00451
00452 void UnitTests::testFlashData()
00453 {
00454 FlashData flash(2,3);
00455 double dd[2][3]={{5, 4, 3},
00456 {8, 10, 20}};
00457
00458
00459 double (*pd)[3];
00460 double nT;
00461 VecDouble vComp(3);
00462 VecDouble phases(2);
00463 Matrix phasesC(2,3);
00464 pd=dd;
00465 flash.setData(&(dd[0][0]));
00466 TEST(flash.getMoles(0,0)==5);
00467 TEST(flash.getMoles(0,1)==4);
00468 TEST(flash.getMoles(0,2)==3);
00469 TEST(flash.getMoles(1,0)==8);
00470 TEST(flash.getMoles(1,1)==10);
00471 TEST(flash.getMoles(1,2)==20);
00472 TEST(flash.getTotalMoles()==50.0);
00473 flash.getMixtureComposition(vComp,nT);
00474 TEST(nT==50);
00475 TEST(NumericMethods::a_equal(vComp(0),13.0/50.0,1.e-12));
00476 TEST(NumericMethods::a_equal(vComp(1),14.0/50.0,1.e-12));
00477 TEST(NumericMethods::a_equal(vComp(2),23.0/50.0,1.e-12));
00478 flash.getPhaseConcentrationsAndTotalMoles(phasesC,phases);
00479 TEST(NumericMethods::a_equal(phasesC(0,0),5.0/12.0,1.e-12));
00480 TEST(NumericMethods::a_equal(phasesC(0,1),4.0/12.0,1.e-12));
00481 TEST(NumericMethods::a_equal(phasesC(0,2),3.0/12.0,1.e-12));
00482 TEST(NumericMethods::a_equal(phasesC(1,0),8.0/38.0,1.e-12));
00483 TEST(NumericMethods::a_equal(phasesC(1,1),10.0/38.0,1.e-12));
00484 TEST(NumericMethods::a_equal(phasesC(1,2),20.0/38.0,1.e-12));
00485 TEST(NumericMethods::a_equal(phases(0),12.0,1.e-12));
00486 TEST(NumericMethods::a_equal(phases(1),38.0,1.e-12));
00487
00488
00489 VecDouble vTC(3);
00490
00491 flash.getTotalComponentsMoles(vTC);
00492 TEST(vTC(0)==13);
00493 TEST(vTC(1)==14);
00494 TEST(vTC(2)==23);
00495
00496
00497 VecDouble X(3);
00498 VecDouble resp(2);
00499 X(0)=1;
00500 X(1)=-2;
00501 X(2)=3;
00502 flash.multiply(resp,X);
00503 TEST(NumericMethods::a_equal(resp(0),6,1.e-12));
00504 TEST(NumericMethods::a_equal(resp(1),48,1.e-12));
00505
00506
00507 flash.zeroEntries();
00508 TEST(flash.getMoles(0,0)==0);
00509 TEST(flash.getMoles(0,1)==0);
00510 TEST(flash.getMoles(0,2)==0);
00511 TEST(flash.getMoles(1,0)==0);
00512 TEST(flash.getMoles(1,1)==0);
00513 TEST(flash.getMoles(1,2)==0);
00514
00515
00516
00517
00518
00519 }
00520
00521
00522 void UnitTests::testFlashCO2Brine()
00523 {
00524 int SALT= FlashCO2Brine::SALT;
00525 int WATER= FlashCO2Brine::WATER;
00526 int CO2= FlashCO2Brine::CO2;
00527 int AQUEOUS= FlashCO2Brine::AQUEOUS;
00528 int CO2_RICH= FlashCO2Brine::CO2_RICH;
00529
00530 Flash::BrineCO2Data data;
00531 Flash::GeometricFlashBrineCO2(data,353,270,0.1,1,1);
00532
00533 OrthoMesh mesh(Point3D(0,0,0),Point3D(1,1,1),10,10,10);
00534 FILE *fd;
00535
00536
00537 fd=fopen("Figure5-1.dat","wr");
00538
00539 FlashData data1(2,3);
00540 VecDouble moles(3);
00541
00542 FlashCO2Brine flash20(mesh,Units::CtoK(20));
00543 FlashCO2Brine flash25(mesh,Units::CtoK(25));
00544 FlashCO2Brine flash(mesh,Units::CtoK(30));
00545 FlashCO2Brine flash90(mesh,Units::CtoK(90));
00546 FlashCO2Brine flash40(mesh,Units::CtoK(40));
00547 FlashCO2Brine flash60(mesh,Units::CtoK(60));
00548 FlashCO2Brine flash80(mesh,Units::CtoK(80));
00549 double mError=0.0;
00550 moles(2)=0;
00551 printf("Volume %g %g\n",ComputeVolumeCO2(303,72,Flash::GAS),ComputeVolumeCO2(303,72,Flash::LIQUID));
00552 for (double P=60e+5;P<600e+5;P+=10e+5)
00553 {
00554 for (double no=0.0;no<100;no++)
00555 {
00556 printf("Volume %g %g\n",ComputeVolumeCO2(Units::CtoK(25),Units::NewtonToBars(P),Flash::GAS),ComputeVolumeCO2(Units::CtoK(25),Units::NewtonToBars(P),Flash::LIQUID));
00557 double ng=100-no;
00558 moles(0)=no;
00559 moles(1)=ng;
00560 flash25.flash(P,moles,data1);
00561 data1.print();
00562 double b = flash25.getFluidCompressibility(P,data1);
00563 double c = flash25.getNumericFluidCompressibility(P,data1,1.e-7);
00564 printf("Error P: %6g <> %6g = %6g\n",b,c,fabs((b-c)/b));
00565 if (mError < (fabs(b-c)/b ))
00566 {
00567 mError = fabs(b-c)/b;
00568 }
00569 }
00570
00571 }
00572 printf("Max Error: %g\n",mError);
00573
00574
00575
00576 moles(CO2)=10;
00577 moles(WATER)=20;
00578 fprintf(fd,"#Figure 5.1 for macl = 0,1,2,3,4\n");
00579 for (double P=1e+5;P<600e+5;P+=10e+5)
00580 {
00581 fprintf(fd,"%g ",P);
00582 for (unsigned i=0;i<5;i++)
00583 {
00584 moles(SALT)=i*moles(WATER)*flash.getComponentsMolarMass()(WATER);
00585 flash.flash(P,moles,data1);
00586 fprintf(fd," %g",data1.getMoles(AQUEOUS,CO2)/(moles(WATER)*flash.getComponentsMolarMass()(WATER)));
00587 }
00588 fprintf(fd,"\n");
00589
00590
00591 }
00592 fclose(fd);
00593
00594 fd=fopen("Figure5-3.dat","wr");
00595 fprintf(fd,"#Figure 5.3 for macl = 0,1,2,3,4\n");
00596 for (double P=1e+5;P<=600e+5;P+=10e+5)
00597 {
00598 fprintf(fd,"%g ",P);
00599 for (unsigned i=0;i<5;i++)
00600 {
00601 moles(SALT)=i*moles(WATER)*flash90.getComponentsMolarMass()(WATER);
00602 flash90.flash(P,moles,data1);
00603 fprintf(fd," %g",data1.getMoles(AQUEOUS,CO2)/(moles(WATER)*flash90.getComponentsMolarMass()(WATER)));
00604 }
00605 fprintf(fd,"\n");
00606
00607
00608
00609 }
00610 fclose(fd);
00611
00612
00613 fd=fopen("Figure5-4.dat","wr");
00614 fprintf(fd,"#Figure 5.4 Yh2ofor macl = 0,1,2,3,4\n");
00615 for (double P=1e+5;P<200e+5;P+=1e+5)
00616 {
00617 fprintf(fd,"%g ",P);
00618 for (unsigned i=0;i<5;i++)
00619 {
00620 moles(SALT)=i*moles(WATER)*flash.getComponentsMolarMass()(WATER);
00621 flash.flash(P,moles,data1);
00622 fprintf(fd," %g",1000*data1.getMoles(CO2_RICH,WATER)/(data1.getMoles(CO2_RICH,WATER)+ data1.getMoles(CO2_RICH,CO2)));
00623 }
00624 fprintf(fd,"\n");
00625 }
00626 fclose(fd);
00627
00628 fd=fopen("Figure5-6.dat","wr");
00629 fprintf(fd,"#Figure 5.6 Yh2ofor macl = 0,1,2,3,4\n");
00630 for (double P=1e+5;P<200e+5;P+=1e+5)
00631 {
00632 fprintf(fd,"%g ",P);
00633 for (unsigned i=0;i<5;i++)
00634 {
00635 moles(SALT)=i*moles(WATER)*flash90.getComponentsMolarMass()(WATER);
00636 flash90.flash(P,moles,data1);
00637 fprintf(fd," %g",1000*data1.getMoles(CO2_RICH,WATER)/(data1.getMoles(CO2_RICH,WATER)+ data1.getMoles(CO2_RICH,CO2)));
00638 }
00639 fprintf(fd,"\n");
00640 }
00641 fclose(fd);
00642
00643 fd=fopen("Figure5-7.dat","wr");
00644 fprintf(fd,"#Figure 5.7 mco2/KgH20 T = 40,60,80\n");
00645 for (double P=1e+5;P<100e+5;P+=1e+5)
00646 {
00647 fprintf(fd,"%g ",P);
00648
00649 moles(SALT)=3.997*moles(WATER)*flash40.getComponentsMolarMass()(WATER);
00650 flash40.flash(P,moles,data1);
00651 fprintf(fd," %g",data1.getMoles(AQUEOUS,CO2)/(data1.getMoles(AQUEOUS,WATER)*flash40.getComponentsMolarMass()(WATER)));
00652
00653 moles(SALT)=3.997*moles(WATER)*flash60.getComponentsMolarMass()(WATER);
00654 flash60.flash(P,moles,data1);
00655 fprintf(fd," %g",data1.getMoles(AQUEOUS,CO2)/(data1.getMoles(AQUEOUS,WATER)*flash60.getComponentsMolarMass()(WATER)));
00656
00657 moles(SALT)=4.001*moles(WATER)*flash80.getComponentsMolarMass()(WATER);
00658 flash80.flash(P,moles,data1);
00659 fprintf(fd," %g",data1.getMoles(AQUEOUS,CO2)/(data1.getMoles(AQUEOUS,WATER)*flash80.getComponentsMolarMass()(WATER)));
00660
00661
00662 fprintf(fd,"\n");
00663 }
00664 fclose(fd);
00665
00666
00667
00668 fd=fopen("Figure5-12.dat","wr");
00669 fprintf(fd,"#Figure 5.12 mco2/KgH20 T = 40,60,80\n");
00670 moles(WATER)=60;
00671 moles(CO2)=35;
00672 moles(SALT)=5;
00673 for (double P=1e+5;P<100e+5;P+=1e+5)
00674 {
00675 fprintf(fd,"%g ",P);
00676
00677 flash20.flash(P,moles,data1);
00678 fprintf(fd," %g",data1.getMoles(AQUEOUS,CO2)/(data1.getMoles(AQUEOUS,WATER)*flash40.getComponentsMolarMass()(WATER)));
00679
00680
00681 flash25.flash(P,moles,data1);
00682 fprintf(fd," %g",data1.getMoles(AQUEOUS,CO2)/(data1.getMoles(AQUEOUS,WATER)*flash40.getComponentsMolarMass()(WATER)));
00683
00684 flash.flash(P,moles,data1);
00685 fprintf(fd," %g",data1.getMoles(AQUEOUS,CO2)/(data1.getMoles(AQUEOUS,WATER)*flash40.getComponentsMolarMass()(WATER)));
00686
00687 flash40.flash(P,moles,data1);
00688 fprintf(fd," %g",data1.getMoles(AQUEOUS,CO2)/(data1.getMoles(AQUEOUS,WATER)*flash40.getComponentsMolarMass()(WATER)));
00689
00690 flash60.flash(P,moles,data1);
00691 fprintf(fd," %g",data1.getMoles(AQUEOUS,CO2)/(data1.getMoles(AQUEOUS,WATER)*flash40.getComponentsMolarMass()(WATER)));
00692
00693 flash80.flash(P,moles,data1);
00694 fprintf(fd," %g",data1.getMoles(AQUEOUS,CO2)/(data1.getMoles(AQUEOUS,WATER)*flash40.getComponentsMolarMass()(WATER)));
00695
00696 fprintf(fd,"\n");
00697 }
00698 fclose(fd);
00699
00700
00701 fd=fopen("Figure5-13.dat","wr");
00702 fprintf(fd,"#Figure 5.13 mco2/KgH20 T = 40,60,80\n");
00703 moles(WATER)=60;
00704 moles(CO2)=35;
00705 moles(SALT)=5;
00706 double P=Units::BarsToNewton(80);
00707 for (moles(CO2)=0;moles(CO2)<=100;moles(CO2)+=1)
00708 {
00709 fprintf(fd,"%g ",moles(CO2));
00710
00711 moles(SALT)=0;
00712 moles(WATER)=100-moles(CO2)-moles(SALT);
00713 if (moles(WATER)<=0.0)
00714 {
00715 moles(SALT)=100-moles(CO2);
00716 moles(WATER)=0.0;
00717 }
00718 flash25.flash(P,moles,data1);
00719 if (data1.getMoles(AQUEOUS,CO2) == 0.0)
00720 fprintf(fd," 0");
00721 else
00722 fprintf(fd," %g",data1.getMoles(AQUEOUS,CO2)/(data1.getMoles(AQUEOUS,WATER)*flash40.getComponentsMolarMass()(WATER)));
00723
00724 moles(SALT)=0.4;
00725 moles(WATER)=100-moles(CO2)-moles(SALT);
00726 if (moles(WATER)<=0.0)
00727 {
00728 moles(SALT)=100-moles(CO2);
00729 moles(WATER)=0.0;
00730 }
00731 flash25.flash(P,moles,data1);
00732 if (data1.getMoles(AQUEOUS,CO2) == 0.0)
00733 fprintf(fd," 0");
00734 else
00735 fprintf(fd," %g",data1.getMoles(AQUEOUS,CO2)/(data1.getMoles(AQUEOUS,WATER)*flash40.getComponentsMolarMass()(WATER)));
00736
00737
00738 moles(SALT)=0.8;
00739 moles(WATER)=100-moles(CO2)-moles(SALT);
00740 if (moles(WATER)<=0.0)
00741 {
00742 moles(SALT)=100-moles(CO2);
00743 moles(WATER)=0.0;
00744 }
00745 flash25.flash(P,moles,data1);
00746 if (data1.getMoles(AQUEOUS,CO2) == 0.0)
00747 fprintf(fd," 0");
00748 else
00749 fprintf(fd," %g",data1.getMoles(AQUEOUS,CO2)/(data1.getMoles(AQUEOUS,WATER)*flash40.getComponentsMolarMass()(WATER)));
00750
00751
00752 moles(SALT)=1.2;
00753 moles(WATER)=100-moles(CO2)-moles(SALT);
00754 if (moles(WATER)<=0.0)
00755 {
00756 moles(SALT)=100-moles(CO2);
00757 moles(WATER)=0.0;
00758 }
00759 flash25.flash(P,moles,data1);
00760 if (data1.getMoles(AQUEOUS,CO2) == 0.0)
00761 fprintf(fd," 0");
00762 else
00763 fprintf(fd," %g",data1.getMoles(AQUEOUS,CO2)/(data1.getMoles(AQUEOUS,WATER)*flash40.getComponentsMolarMass()(WATER)));
00764
00765 moles(SALT)=1.6;
00766 moles(WATER)=100-moles(CO2)-moles(SALT);
00767 if (moles(WATER)<=0.0)
00768 {
00769 moles(SALT)=100-moles(CO2);
00770 moles(WATER)=0.0;
00771 }
00772 flash25.flash(P,moles,data1);
00773 if (data1.getMoles(AQUEOUS,CO2) == 0.0)
00774 fprintf(fd," 0");
00775 else
00776 fprintf(fd," %g",data1.getMoles(AQUEOUS,CO2)/(data1.getMoles(AQUEOUS,WATER)*flash40.getComponentsMolarMass()(WATER)));
00777
00778 moles(SALT)=2.0;
00779 moles(WATER)=100-moles(CO2)-moles(SALT);
00780 if (moles(WATER)<=0.0)
00781 {
00782 moles(SALT)=100-moles(CO2);
00783 moles(WATER)=0.0;
00784 }
00785 flash25.flash(P,moles,data1);
00786 if (data1.getMoles(AQUEOUS,CO2) == 0.0)
00787 fprintf(fd," 0");
00788 else
00789 fprintf(fd," %g",data1.getMoles(AQUEOUS,CO2)/(data1.getMoles(AQUEOUS,WATER)*flash40.getComponentsMolarMass()(WATER)));
00790
00791
00792 fprintf(fd,"\n");
00793 }
00794 fclose(fd);
00795
00796
00797
00798 fd=fopen("Figure5-14.dat","wr");
00799 fprintf(fd,"#Figure 5.14 mco2/KgH20 T = 40,60,80\n");
00800 moles(WATER)=60;
00801 moles(CO2)=35;
00802 moles(SALT)=5;
00803 P=Units::BarsToNewton(80);
00804 for (moles(CO2)=0;moles(CO2)<=100;moles(CO2)+=1)
00805 {
00806 fprintf(fd,"%g ",moles(CO2));
00807
00808 moles(SALT)=0;
00809 moles(WATER)=100-moles(CO2)-moles(SALT);
00810 if (moles(WATER)<=0.0)
00811 {
00812 moles(SALT)=100-moles(CO2);
00813 moles(WATER)=0.0;
00814 }
00815 flash60.flash(P,moles,data1);
00816 if (data1.getMoles(AQUEOUS,CO2) == 0.0)
00817 fprintf(fd," 0");
00818 else
00819 fprintf(fd," %g",data1.getMoles(AQUEOUS,CO2)/(data1.getMoles(AQUEOUS,WATER)*flash40.getComponentsMolarMass()(WATER)));
00820
00821 moles(SALT)=0.4;
00822 moles(WATER)=100-moles(CO2)-moles(SALT);
00823 if (moles(WATER)<=0.0)
00824 {
00825 moles(SALT)=100-moles(CO2);
00826 moles(WATER)=0.0;
00827 }
00828 flash60.flash(P,moles,data1);
00829 if (data1.getMoles(AQUEOUS,CO2) == 0.0)
00830 fprintf(fd," 0");
00831 else
00832 fprintf(fd," %g",data1.getMoles(AQUEOUS,CO2)/(data1.getMoles(AQUEOUS,WATER)*flash40.getComponentsMolarMass()(WATER)));
00833
00834
00835 moles(SALT)=0.8;
00836 moles(WATER)=100-moles(CO2)-moles(SALT);
00837 if (moles(WATER)<=0.0)
00838 {
00839 moles(SALT)=100-moles(CO2);
00840 moles(WATER)=0.0;
00841 }
00842 flash60.flash(P,moles,data1);
00843 if (data1.getMoles(AQUEOUS,CO2) == 0.0)
00844 fprintf(fd," 0");
00845 else
00846 fprintf(fd," %g",data1.getMoles(AQUEOUS,CO2)/(data1.getMoles(AQUEOUS,WATER)*flash40.getComponentsMolarMass()(WATER)));
00847
00848
00849 moles(SALT)=1.2;
00850 moles(WATER)=100-moles(CO2)-moles(SALT);
00851 if (moles(WATER)<=0.0)
00852 {
00853 moles(SALT)=100-moles(CO2);
00854 moles(WATER)=0.0;
00855 }
00856 flash60.flash(P,moles,data1);
00857 if (data1.getMoles(AQUEOUS,CO2) == 0.0)
00858 fprintf(fd," 0");
00859 else
00860 fprintf(fd," %g",data1.getMoles(AQUEOUS,CO2)/(data1.getMoles(AQUEOUS,WATER)*flash40.getComponentsMolarMass()(WATER)));
00861
00862 moles(SALT)=1.6;
00863 moles(WATER)=100-moles(CO2)-moles(SALT);
00864 if (moles(WATER)<=0.0)
00865 {
00866 moles(SALT)=100-moles(CO2);
00867 moles(WATER)=0.0;
00868 }
00869 flash60.flash(P,moles,data1);
00870 if (data1.getMoles(AQUEOUS,CO2) == 0.0)
00871 fprintf(fd," 0");
00872 else
00873 fprintf(fd," %g",data1.getMoles(AQUEOUS,CO2)/(data1.getMoles(AQUEOUS,WATER)*flash40.getComponentsMolarMass()(WATER)));
00874
00875 moles(SALT)=2.0;
00876 moles(WATER)=100-moles(CO2)-moles(SALT);
00877 if (moles(WATER)<=0.0)
00878 {
00879 moles(SALT)=100-moles(CO2);
00880 moles(WATER)=0.0;
00881 }
00882 flash60.flash(P,moles,data1);
00883 if (data1.getMoles(AQUEOUS,CO2) == 0.0)
00884 fprintf(fd," 0");
00885 else
00886 fprintf(fd," %g",data1.getMoles(AQUEOUS,CO2)/(data1.getMoles(AQUEOUS,WATER)*flash40.getComponentsMolarMass()(WATER)));
00887
00888
00889 fprintf(fd,"\n");
00890 }
00891 fclose(fd);
00892
00893
00894 fd=fopen("Figure5-15.dat","wr");
00895 fprintf(fd,"#Figure 5.15 mco2/KgH20 T = 40,60,80\n");
00896 moles(WATER)=60;
00897 moles(CO2)=35;
00898 moles(SALT)=5;
00899 P=Units::BarsToNewton(80);
00900 double znacl[6]={0,0.4,0.8,1.2,1.6,2.0};
00901 for (moles(CO2)=0;moles(CO2)<=100;moles(CO2)+=1)
00902 {
00903 fprintf(fd,"%g ",moles(CO2));
00904 for (unsigned i=0;i<6;i++)
00905 {
00906 moles(SALT)=znacl[i];
00907 moles(WATER)=100-moles(CO2)-moles(SALT);
00908 if (moles(WATER)<=0.0)
00909 {
00910 moles(SALT)=100-moles(CO2);
00911 moles(WATER)=0.0;
00912 }
00913 flash25.flash(P,moles,data1);
00914 if (data1.getMoles(CO2_RICH,WATER) == 0.0)
00915 fprintf(fd," 0");
00916 else
00917 fprintf(fd," %g",1000*data1.getMoles(CO2_RICH,WATER)/(data1.getMoles(CO2_RICH,WATER)+data1.getMoles(CO2_RICH,CO2)));
00918 }
00919 fprintf(fd,"\n");
00920 }
00921 fclose(fd);
00922
00923
00924 fd=fopen("Figure5-16.dat","wr");
00925 fprintf(fd,"#Figure 5.16 mco2/KgH20 T = 40,60,80\n");
00926 moles(WATER)=60;
00927 moles(CO2)=35;
00928 moles(SALT)=5;
00929 P=Units::BarsToNewton(80);
00930 for (moles(CO2)=0;moles(CO2)<=100;moles(CO2)+=1)
00931 {
00932 fprintf(fd,"%g ",moles(CO2));
00933 for (unsigned i=0;i<6;i++)
00934 {
00935 moles(SALT)=znacl[i];
00936 moles(WATER)=100-moles(CO2)-moles(SALT);
00937 if (moles(WATER)<=0.0)
00938 {
00939 moles(SALT)=100-moles(CO2);
00940 moles(WATER)=0.0;
00941 }
00942 flash60.flash(P,moles,data1);
00943 if (data1.getMoles(CO2_RICH,WATER) == 0.0)
00944 fprintf(fd," 0");
00945 else
00946 fprintf(fd," %g",1000*data1.getMoles(CO2_RICH,WATER)/(data1.getMoles(CO2_RICH,WATER)+data1.getMoles(CO2_RICH,CO2)));
00947 }
00948 fprintf(fd,"\n");
00949 }
00950 fclose(fd);
00951
00952
00953 VecDouble phasesVol(2);
00954 fd = fopen("CO2VolumeCurve.dat","wr");
00955 fprintf(fd,"#Volume Curve as function of pressure (BARS) T = 40,60,80\n");
00956 moles(WATER)=0;
00957 moles(CO2)=1;
00958 moles(SALT)=0;
00959 P=Units::BarsToNewton(80);
00960 data1.zeroEntries();
00961 data1.setMoles(CO2_RICH,CO2,1);
00962 for (P =10;P <= 600; P+=10)
00963 {
00964 fprintf(fd,"%g\t",P);
00965
00966 flash40.getPhasesVolume(Units::BarsToNewton(P),data1,phasesVol);
00967 fprintf(fd,"%g\t",phasesVol(1));
00968
00969 flash60.getPhasesVolume(Units::BarsToNewton(P),data1,phasesVol);
00970 fprintf(fd,"%g\t",phasesVol(1));
00971
00972 flash80.getPhasesVolume(Units::BarsToNewton(P),data1,phasesVol);
00973 fprintf(fd,"%g\t\n",phasesVol(1));
00974 }
00975 fclose(fd);
00976
00977
00978
00979
00980
00981
00982
00983
00984
00985
00986
00987 }
00988
00989
00990
00991
00992 void UnitTests::testSimpleBlackOilFlash()
00993 {
00994 VecDouble v(2),dVt(2),dVtNum(2);
00995 OrthoMesh &mesh = app.getMesh();
00996 FlashSimpleBlackOil flash(mesh);
00997
00998 FlashData data(2,2);
00999 app.printBeginSection("Testing Simple Black Oil Flash");
01000 std::ofstream file("SimpleBlackOil.dat");
01001 file << "#Derivatives Check data\n";
01002 file << "#P <no> <Exact dVt/dP> <Exact dVt/dno> <Exact Dvt/dng> <Approx dVt/dP> <Approx dVt/dno> <Approx Dvt/dng>\n";
01003
01004
01005
01006 for (double P=1; P <1500; P+=10)
01007 {
01008 for (double no=0.0;no<100;no++)
01009 {
01010 file << P << "\t";
01011 file << no << "\t";
01012 double ng=100-no;
01013 v(0)=no;
01014 v(1)=ng;
01015 flash.flash(P,v,data);
01016 data.print();
01017 double b = flash.getFluidCompressibility(P,data);
01018 file << b << "\t";
01019 double c = flash.getNumericFluidCompressibility(P,data,1.e-6);
01020 printf("Error P: %6g <> %6g = %6g\n",b,c,fabs((b-c)/b));
01021
01022 flash.getTotalVolumeDerivatives(P,data,dVt);
01023 flash.getNumericTotalVolumeDerivative(P,v,dVtNum,1.e-5);
01024 if (data.getMoles(FlashSimpleBlackOil::GAS,FlashSimpleBlackOil::GAS) == 0.0)
01025 printf("UNSATURATED\n");
01026 else
01027 printf("SATURATED\n");
01028 printf("Error no: %6g <> %6g = %6g\n",dVt(0),dVtNum(0),fabs((dVt(0)-dVtNum(0))/dVt(0)));
01029 printf("Error ng: %6g <> %6g = %6g\n",dVt(1),dVtNum(1),fabs((dVt(1)-dVtNum(1))/dVt(1)));
01030 file << dVt(0) << "\t" << dVt(1) << "\t" << c << "\t" << dVtNum(0) << "\t" << dVtNum(1) ;
01031 file << std::endl;
01032 }
01033 file << "\n";
01034 }
01035 file.close();
01036
01037 }
01038
01039
01040 void UnitTests::plotMobilities()
01041 {
01042 GeneralFunctionInterface *fMob,*fMobT,*fFracFlux,*fGrav,*DfFracFlux,*DfGrav;
01043 app.getMobilities(fMob,fMobT,fFracFlux,DfFracFlux,fGrav,DfGrav);
01044
01045 GeneralFunctionInterface *fKr,*fDKr;
01046
01047 VecDouble p1(1);
01048 VecDouble p2(1);
01049 p1(0)=0;
01050 p2(0)=1;
01051
01052
01053
01054 std::cout << "Printing Fractional Flow into file fracFlow.dat: ";
01055 Function1D *fFrac1D = dynamic_cast<Function1D*>( fFracFlux);
01056 std::ofstream file("fracFlow.dat");
01057 fFrac1D->plot(file,0,1,2000);
01058 file.close();
01059 std::cout << "done" << std::endl;
01060
01061 file.open("fracFlowD.dat");
01062 std::cout << "Printing Fractional Flow Derivative into file fracFlowD.dat: ";
01063 DfFracFlux->plotInLine(p1,p2,2000,file);
01064 std::cout << "done" << std::endl;
01065 file.close();
01066
01067
01068 file.open("TotalMobility.dat");
01069 std::cout << "Printing Total Mobility into file TotalMobility.dat: ";
01070 fMobT->plotInLine(p1,p2,2000,file);
01071 std::cout << "done" << std::endl;
01072 file.close();
01073
01074
01075 file.open("GravMobility.dat");
01076 Function1D *pfGrav = dynamic_cast<Function1D*>( fGrav);
01077 std::cout << "Printing Gravity Mobility into file GravMobility.dat: ";
01078 pfGrav->plot(file,0,1,2000);
01079 std::cout << "done" << std::endl;
01080 file.close();
01081
01082 file.open("DGravMobility.dat");
01083 std::cout << "Printing Derivative Gravity Mobility into file DGravMobility.dat: ";
01084 Function1D *pDfGrav = dynamic_cast<Function1D*>( DfGrav);
01085 pDfGrav->plot(file,0,1,2000);
01086 std::cout << "done" << std::endl;
01087 file.close();
01088
01089 app.getKr(fKr,fDKr);
01090 if (fKr)
01091 {
01092 for (unsigned cmp=0;cmp < fKr->getImageDim();cmp++)
01093 {
01094 char fileName[100];
01095 sprintf(fileName,"fKr_%u.dat",cmp);
01096 file.open(fileName);
01097 std::cout << "Printing Kr" << cmp << " into file " << fileName << ".....";
01098 fKr->plotInLine(p1,p2,2000,file,cmp);
01099 std::cout << "Ok" << std::endl;
01100 file.close();
01101 }
01102 }
01103 if (fDKr)
01104 {
01105 for (unsigned cmp=0;cmp < fDKr->getImageDim();cmp++)
01106 {
01107 char fileName[100];
01108 sprintf(fileName,"fDKr_%u.dat",cmp);
01109 file.open(fileName);
01110 std::cout << "Printing Dx(Kr" << cmp << ") into file " << fileName << ".....";
01111 fDKr->plotInLine(p1,p2,2000,file,cmp);
01112 std::cout << "Ok" << std::endl;
01113 file.close();
01114 }
01115 }
01116
01117
01118
01119
01120 }
01121
01122
01123 void UnitTests::plotPc()
01124 {
01125 GeneralFunctionInterface *fPc,*dfPc,*fPcInv;
01126 app.getPcFunction(&fPc,&dfPc,&fPcInv);
01127 std::ofstream file;
01128 Function1D& pfPc = dynamic_cast<Function1D&>(*fPc);
01129 Function1D& pdfPc = dynamic_cast<Function1D&>(*dfPc);
01130
01131
01132
01133 file.open("Pc.dat");
01134 std::cout << "Printing Cappillary Pressure curve to Pc.dat ..................";
01135 pfPc.plot(file,0,1,2000);
01136 std::cout << "done" << std::endl;
01137 file.close();
01138
01139 file.open("DPc.dat");
01140 std::cout << "Printing Derivative of Cappillary Pressure curve to DPc.dat ....";
01141 pdfPc.plot(file,0,1,2000);
01142 std::cout << "done" << std::endl;
01143 file.close();
01144
01145 if (fPcInv)
01146 {
01147 Function1D* pfPcInv = dynamic_cast<Function1D*>(fPcInv);
01148 if (pfPcInv)
01149 {
01150 file.open("PcInv.dat");
01151 std::cout << "Printing Inverse of Cappillary Pressure at PcInv.dat ....";
01152 pfPcInv->plot(file,0,1,2000);
01153 std::cout << "done" << std::endl;
01154 file.close();
01155 }
01156 }
01157
01158
01159
01160 }
01161
01162 void getValue(std::string str,double *d)
01163 {
01164 std::cout << str.c_str() << " ";
01165 std::cin >> *d;
01166
01167 }
01168
01169
01170 void UnitTests::calculateFlashCO2Brine()
01171 {
01172 VecDouble v(3);
01173 double T,P;
01174 printf("Type the components compositions:");
01175 getValue("Water:", &v(0));
01176 getValue("CO2:", &v(1));
01177 getValue("Salt:", &v(2));
01178 getValue("Temperature: ",&T);
01179 getValue("Pressure: ",&P);
01180
01181 OrthoMesh &mesh = app.getMesh();
01182 FlashCO2Brine flash(mesh,T);
01183 VecDouble phasesVol(2),phasesDens(2),phasesT(2);
01184 FlashData data(2,3);
01185 flash.flash(P,v,data);
01186 flash.getPhasesVolume(P,data,phasesVol);
01187 data.getPhasesDensities(phasesDens,flash.getComponentsMolarMass(),phasesVol);
01188 data.getPhasesTotalMoles(phasesT);
01189 std::cout << "Solution:\n";
01190 std::cout.precision(10);
01191 data.print();
01192 std::cout.setf(std::ios::scientific,std::ios::floatfield);
01193 std::cout << "Volumes: " << phasesVol << std::endl;
01194 std::cout << "Densities: " << phasesDens << std::endl;
01195 std::cout << "CO2 supercritic molar volume: " << flash.getPureCO2MolarVolume(P,T) << std::endl;
01196 std::cout << "CO2 supercritic density (kg/m3): " <<flash.getComponentsMolarMass()(FlashCO2Brine::CO2)/ flash.getPureCO2MolarVolume(P,T) << std::endl;
01197 std::cout << "Total Moles in each phase: " << phasesT;
01198 std::cout << std::endl;
01199
01200 }
