diff --git a/benchmarks/upscale_relperm_benchmark.cpp b/benchmarks/upscale_relperm_benchmark.cpp index 9cd9aa3..a0d4964 100644 --- a/benchmarks/upscale_relperm_benchmark.cpp +++ b/benchmarks/upscale_relperm_benchmark.cpp @@ -189,7 +189,6 @@ try // Variables used for timing/profiling: clock_t start, finish; double timeused = 0.0, timeused_tesselation = 0.0; - double timeused_upscale_wallclock = 0.0; clock_t global_start = clock(); // Timing used for benchmarking @@ -295,11 +294,9 @@ try finish = clock(); timeused = (double(finish)-double(start))/CLOCKS_PER_SEC; if (helper.isMaster) cout << " (" << timeused <<" secs)" << endl; - const double maxPermContrast = atof(options["maxPermContrast"].c_str()); const double minPerm = atof(options["minPerm"].c_str()); const double maxPerm = atof(options["maxPerm"].c_str()); const double minPoro = atof(options["minPoro"].c_str()); - const double saturationThreshold = atof(options["saturationThreshold"].c_str()); start = clock(); helper.sanityCheckInput(deck, minPerm, maxPerm, minPoro); @@ -454,9 +451,6 @@ try helper.upscaleCapillaryPressure(options, dP); - MonotCubicInterpolator CapPressureVsWaterSaturation(helper.WaterSaturationVsCapPressure.get_fVector(), - helper.WaterSaturationVsCapPressure.get_xVector()); - clock_t start_upscaling = clock(); /***************************************************************************** @@ -469,7 +463,6 @@ try */ helper.upscaleSinglePhasePermeability(); - typedef SinglePhaseUpscaler::permtensor_t Matrix; /***************************************************************** * Step 8: @@ -487,262 +480,9 @@ try * c: Calculate relperm tensors from all the phase perm tensors. */ - // Put correct number of zeros in, just to be able to access RelPerm[index] later - helper.WaterSaturation.resize(helper.points, 0.0); - for (size_t i = 0; i < (helper.upscaleBothPhases?2:1); ++i) - helper.PhasePerm[i].resize(helper.points, std::vector(helper.tensorElementCount)); - - // Make vector of capillary pressure points corresponding to uniformly distribued - // saturation points between Swor and Swir. - - for (int pointidx = 1; pointidx <= helper.points; ++pointidx) { - // pointidx=1 corresponds to Swir, pointidx=points to Swor. - double saturation = helper.Swir + (helper.Swor-helper.Swir)/(helper.points-1)*(pointidx-1); - helper.pressurePoints.push_back(CapPressureVsWaterSaturation.evaluate(saturation)); - } - // Preserve max and min pressures - helper.pressurePoints[0]=helper.Pcmax; - helper.pressurePoints[helper.pressurePoints.size()-1]=helper.Pcmin; - - // Fill with zeros initially (in case of non-mpi) - helper.node_vs_pressurepoint.resize(helper.points); - -#if HAVE_MPI - // Distribute work load over mpi nodes. - for (int idx=0; idx < points; ++idx) { - // Ensure master node gets equal or less work than the other nodes, since - // master node also computes single phase perm. - node_vs_pressurepoint[idx] = (mpi_nodecount-1) - idx % mpi_nodecount; - /*if (helper.isMaster) { - cout << "Pressure point " << idx << " assigned to node " << node_vs_pressurepoint[idx] << endl; - }*/ - } -#endif - - - clock_t start_upscale_wallclock = clock(); - - double waterVolumeLF; - // Now loop through the vector of capillary pressure points that - // this node should compute. - for (int pointidx = 0; pointidx < helper.points; ++pointidx) { - - // Should "I" (mpi-wise) compute this pressure point? - if (helper.node_vs_pressurepoint[pointidx] == mpi_rank) { - - double Ptestvalue = helper.pressurePoints[pointidx]; - - double accPhasePerm = 0.0; - double accPhase2Perm = 0.0; - - double maxPhasePerm = 0.0; - double maxPhase2Perm = 0.0; - - vector phasePermValues, phase2PermValues; - vector > phasePermValuesDiag, phase2PermValuesDiag; - phasePermValues.resize(helper.satnums.size()); - phasePermValuesDiag.resize(helper.satnums.size()); - if (helper.upscaleBothPhases) { - phase2PermValues.resize(helper.satnums.size()); - phase2PermValuesDiag.resize(helper.satnums.size()); - } - waterVolumeLF = 0.0; - for (unsigned int i = 0; i < ecl_idx.size(); ++i) { - unsigned int cell_idx = ecl_idx[i]; - double cellPhasePerm = minPerm; - double cellPhase2Perm = minPerm; - vector cellPhasePermDiag, cellPhase2PermDiag; - cellPhasePermDiag.push_back(minPerm); - cellPhasePermDiag.push_back(minPerm); - cellPhasePermDiag.push_back(minPerm); - if (helper.upscaleBothPhases) { - cellPhase2PermDiag.push_back(minPerm); - cellPhase2PermDiag.push_back(minPerm); - cellPhase2PermDiag.push_back(minPerm); - } - - if (helper.satnums[cell_idx] > 0) { // handle "no rock" cells with satnum zero - // cout << endl << "Cell no. " << cell_idx << endl; - double PtestvalueCell; - if (includeGravity) { - PtestvalueCell = Ptestvalue - dP[cell_idx]; - } - else { - PtestvalueCell = Ptestvalue; - } - - if (! helper.anisotropic_input) { - - double Jvalue = sqrt(helper.perms[0][cell_idx] * milliDarcyToSqMetre/helper.poros[cell_idx]) * PtestvalueCell; - //cout << "JvalueCell: " << Jvalue << endl; - double WaterSaturationCell - = helper.InvJfunctions[int(helper.satnums[cell_idx])-1].evaluate(Jvalue); - waterVolumeLF += WaterSaturationCell * helper.cellPoreVolumes[cell_idx]; - - // Compute cell relative permeability. We use a lower cutoff-value as we - // easily divide by zero here. When water saturation is - // zero, we get 'inf', which is circumvented by the cutoff value. - cellPhasePerm = - helper.Krfunctions[0][0][int(helper.satnums[cell_idx])-1].evaluate(WaterSaturationCell) * - helper.perms[0][cell_idx]; - if (helper.upscaleBothPhases) { - cellPhase2Perm = - helper.Krfunctions[0][1][int(helper.satnums[cell_idx])-1].evaluate(WaterSaturationCell) * - helper.perms[0][cell_idx]; - } - } - else { - double WaterSaturationCell = helper.SwPcfunctions[int(helper.satnums[cell_idx])-1].evaluate(PtestvalueCell); - //cout << PtestvalueCell << "\t" << helper.WaterSaturationCell << endl; - waterVolumeLF += WaterSaturationCell * helper.cellPoreVolumes[cell_idx]; - - cellPhasePermDiag[0] = helper.Krfunctions[0][0][int(helper.satnums[cell_idx])-1].evaluate(WaterSaturationCell) * - helper.perms[0][cell_idx]; - cellPhasePermDiag[1] = helper.Krfunctions[1][0][int(helper.satnums[cell_idx])-1].evaluate(WaterSaturationCell) * - helper.perms[1][cell_idx]; - cellPhasePermDiag[2] = helper.Krfunctions[2][0][int(helper.satnums[cell_idx])-1].evaluate(WaterSaturationCell) * - helper.perms[2][cell_idx]; - if (helper.upscaleBothPhases) { - cellPhase2PermDiag[0] = helper.Krfunctions[0][1][int(helper.satnums[cell_idx])-1].evaluate(WaterSaturationCell) * - helper.perms[0][cell_idx]; - cellPhase2PermDiag[1] = helper.Krfunctions[1][1][int(helper.satnums[cell_idx])-1].evaluate(WaterSaturationCell) * - helper.perms[1][cell_idx]; - cellPhase2PermDiag[2] = helper.Krfunctions[2][1][int(helper.satnums[cell_idx])-1].evaluate(WaterSaturationCell) * - helper.perms[2][cell_idx]; - } - } - - phasePermValues[cell_idx] = cellPhasePerm; - phasePermValuesDiag[cell_idx] = cellPhasePermDiag; - maxPhasePerm = max(maxPhasePerm, cellPhasePerm); - maxPhasePerm = max(maxPhasePerm, *max_element(cellPhasePermDiag.begin(), - cellPhasePermDiag.end())); - if (helper.upscaleBothPhases) { - phase2PermValues[cell_idx] = cellPhase2Perm; - phase2PermValuesDiag[cell_idx] = cellPhase2PermDiag; - maxPhase2Perm = max(maxPhase2Perm, cellPhase2Perm); - maxPhase2Perm = max(maxPhase2Perm, *max_element(cellPhase2PermDiag.begin(), - cellPhase2PermDiag.end())); - } - } - } - // Now we can determine the smallest permitted permeability we can calculate for - - // We have both a fixed bottom limit, as well as a possible higher limit determined - // by a maximum allowable permeability. - double minPhasePerm = max(maxPhasePerm/maxPermContrast, minPerm); - double minPhase2Perm; - if (helper.upscaleBothPhases) minPhase2Perm = max(maxPhase2Perm/maxPermContrast, minPerm); - - // Now remodel the phase permeabilities obeying minPhasePerm - Matrix cellperm(3,3,nullptr); - zero(cellperm); - for (unsigned int i = 0; i < ecl_idx.size(); ++i) { - unsigned int cell_idx = ecl_idx[i]; - zero(cellperm); - if (! helper.anisotropic_input) { - double cellPhasePerm = max(minPhasePerm, phasePermValues[cell_idx]); - accPhasePerm += cellPhasePerm; - double kval = max(minPhasePerm, cellPhasePerm); - cellperm(0,0) = kval; - cellperm(1,1) = kval; - cellperm(2,2) = kval; - } - else { // anisotropic_input - // Truncate values lower than minPhasePerm upwards. - phasePermValuesDiag[cell_idx][0] = max(minPhasePerm, phasePermValuesDiag[cell_idx][0]); - phasePermValuesDiag[cell_idx][1] = max(minPhasePerm, phasePermValuesDiag[cell_idx][1]); - phasePermValuesDiag[cell_idx][2] = max(minPhasePerm, phasePermValuesDiag[cell_idx][2]); - accPhasePerm += phasePermValuesDiag[cell_idx][0]; // not correct anyway - cellperm(0,0) = phasePermValuesDiag[cell_idx][0]; - cellperm(1,1) = phasePermValuesDiag[cell_idx][1]; - cellperm(2,2) = phasePermValuesDiag[cell_idx][2]; - } - helper.upscaler.setPermeability(i, cellperm); - } - - // Output average phase perm, this is just a reality check so that we are not way off. - //cout << ", Arith. mean phase perm = " << accPhasePerm/float(tesselatedCells) << " mD, "; - - // Call single-phase upscaling code - Matrix phasePermTensor = helper.upscaler.upscaleSinglePhase(); - - // Now upscale phase permeability for phase 2 - Matrix phase2PermTensor; - if (helper.upscaleBothPhases) { - zero(cellperm); - for (unsigned int i = 0; i < ecl_idx.size(); ++i) { - unsigned int cell_idx = ecl_idx[i]; - zero(cellperm); - if (! helper.anisotropic_input) { - double cellPhase2Perm = max(minPhase2Perm, phase2PermValues[cell_idx]); - accPhase2Perm += cellPhase2Perm; - double kval = max(minPhase2Perm, cellPhase2Perm); - cellperm(0,0) = kval; - cellperm(1,1) = kval; - cellperm(2,2) = kval; - } - else { // anisotropic_input - // Truncate values lower than minPhasePerm upwards. - phase2PermValuesDiag[cell_idx][0] = max(minPhase2Perm, phase2PermValuesDiag[cell_idx][0]); - phase2PermValuesDiag[cell_idx][1] = max(minPhase2Perm, phase2PermValuesDiag[cell_idx][1]); - phase2PermValuesDiag[cell_idx][2] = max(minPhase2Perm, phase2PermValuesDiag[cell_idx][2]); - accPhase2Perm += phase2PermValuesDiag[cell_idx][0]; // not correct anyway - cellperm(0,0) = phase2PermValuesDiag[cell_idx][0]; - cellperm(1,1) = phase2PermValuesDiag[cell_idx][1]; - cellperm(2,2) = phase2PermValuesDiag[cell_idx][2]; - } - helper.upscaler.setPermeability(i, cellperm); - } - phase2PermTensor = helper.upscaler.upscaleSinglePhase(); - } - - //cout << phasePermTensor << endl; - - - // Here we recalculate the upscaled water saturation, - // although it is already known when we asked for the - // pressure point to compute for. Nonetheless, we - // recalculate here to avoid any minor roundoff-error and - // interpolation error (this means that the saturation - // points are not perfectly uniformly distributed) - helper.WaterSaturation[pointidx] = waterVolumeLF/helper.poreVolume; - - -#ifdef HAVE_MPI - cout << "Rank " << mpi_rank << ": "; -#endif - cout << Ptestvalue << "\t" << helper.WaterSaturation[pointidx]; - // Store and print phase-perm-result - for (int voigtIdx=0; voigtIdx < helper.tensorElementCount; ++voigtIdx) { - helper.PhasePerm[0][pointidx][voigtIdx] = getVoigtValue(phasePermTensor, voigtIdx); - cout << "\t" << getVoigtValue(phasePermTensor, voigtIdx); - if (helper.upscaleBothPhases){ - helper.PhasePerm[1][pointidx][voigtIdx] = getVoigtValue(phase2PermTensor, voigtIdx); - cout << "\t" << getVoigtValue(phase2PermTensor, voigtIdx); - } - } - cout << endl; - } - } - - clock_t finish_upscale_wallclock = clock(); - timeused_upscale_wallclock = (double(finish_upscale_wallclock)-double(start_upscale_wallclock))/CLOCKS_PER_SEC; - - helper.collectResults(); - - // Average time pr. upscaling point: -#ifdef HAVE_MPI - // Sum the upscaling time used by all processes - double timeused_total; - MPI_Reduce(&timeused_upscale_wallclock, &timeused_total, 1, MPI_DOUBLE, - MPI_SUM, 0, MPI_COMM_WORLD); - double avg_upscaling_time_pr_point = timeused_total/(double)points; - -#else - double avg_upscaling_time_pr_point = timeused_upscale_wallclock / (double)helper.points; -#endif - + double timeused_upscale_wallclock, avg_upscaling_time_pr_point; + std::tie(timeused_upscale_wallclock, avg_upscaling_time_pr_point) = + helper.upscalePermeability(options, dP, mpi_rank); /* * Step 8c: Make relperm values from phaseperms diff --git a/examples/upscale_relperm.cpp b/examples/upscale_relperm.cpp index 3e45d5b..60d6592 100644 --- a/examples/upscale_relperm.cpp +++ b/examples/upscale_relperm.cpp @@ -403,11 +403,9 @@ try res[1] = specgridRecord->getItem("NY")->getInt(0); res[2] = specgridRecord->getItem("NZ")->getInt(0); - const double maxPermContrast = atof(options["maxPermContrast"].c_str()); const double minPerm = atof(options["minPerm"].c_str()); const double maxPerm = atof(options["maxPerm"].c_str()); const double minPoro = atof(options["minPoro"].c_str()); - const double saturationThreshold = atof(options["saturationThreshold"].c_str()); helper.sanityCheckInput(deck, minPerm, maxPerm, minPoro); @@ -673,8 +671,6 @@ try */ helper.upscaleCapillaryPressure(options, dP); - MonotCubicInterpolator CapPressureVsWaterSaturation(helper.WaterSaturationVsCapPressure.get_fVector(), - helper.WaterSaturationVsCapPressure.get_xVector()); /***************************************************************************** * Step 7: @@ -686,7 +682,6 @@ try */ helper.upscaleSinglePhasePermeability(); - typedef SinglePhaseUpscaler::permtensor_t Matrix; /***************************************************************** * Step 8: @@ -704,262 +699,9 @@ try * c: Calculate relperm tensors from all the phase perm tensors. */ - // Put correct number of zeros in, just to be able to access RelPerm[index] later - helper.WaterSaturation.resize(helper.points, 0.0); - for (size_t i = 0; i < (helper.upscaleBothPhases?2:1); ++i) - helper.PhasePerm[i].resize(helper.points, std::vector(helper.tensorElementCount)); - - // Make vector of capillary pressure points corresponding to uniformly distribued - // saturation points between Swor and Swir. - - for (int pointidx = 1; pointidx <= helper.points; ++pointidx) { - // pointidx=1 corresponds to Swir, pointidx=points to Swor. - double saturation = helper.Swir + (helper.Swor-helper.Swir)/(helper.points-1)*(pointidx-1); - helper.pressurePoints.push_back(CapPressureVsWaterSaturation.evaluate(saturation)); - } - // Preserve max and min pressures - helper.pressurePoints[0]=helper.Pcmax; - helper.pressurePoints[helper.pressurePoints.size()-1]=helper.Pcmin; - - // Fill with zeros initially (in case of non-mpi) - helper.node_vs_pressurepoint.resize(helper.points); - -#if HAVE_MPI - // Distribute work load over mpi nodes. - for (int idx=0; idx < points; ++idx) { - // Ensure master node gets equal or less work than the other nodes, since - // master node also computes single phase perm. - node_vs_pressurepoint[idx] = (mpi_nodecount-1) - idx % mpi_nodecount; - /*if (isMaster) { - cout << "Pressure point " << idx << " assigned to node " << node_vs_pressurepoint[idx] << endl; - }*/ - } -#endif - - - clock_t start_upscale_wallclock = clock(); - - double waterVolumeLF; - // Now loop through the vector of capillary pressure points that - // this node should compute. - for (int pointidx = 0; pointidx < helper.points; ++pointidx) { - - // Should "I" (mpi-wise) compute this pressure point? - if (helper.node_vs_pressurepoint[pointidx] == mpi_rank) { - - double Ptestvalue = helper.pressurePoints[pointidx]; - - double accPhasePerm = 0.0; - double accPhase2Perm = 0.0; - - double maxPhasePerm = 0.0; - double maxPhase2Perm = 0.0; - - vector phasePermValues, phase2PermValues; - vector > phasePermValuesDiag, phase2PermValuesDiag; - phasePermValues.resize(helper.satnums.size()); - phasePermValuesDiag.resize(helper.satnums.size()); - if (helper.upscaleBothPhases) { - phase2PermValues.resize(helper.satnums.size()); - phase2PermValuesDiag.resize(helper.satnums.size()); - } - waterVolumeLF = 0.0; - for (unsigned int i = 0; i < ecl_idx.size(); ++i) { - unsigned int cell_idx = ecl_idx[i]; - double cellPhasePerm = minPerm; - double cellPhase2Perm = minPerm; - vector cellPhasePermDiag, cellPhase2PermDiag; - cellPhasePermDiag.push_back(minPerm); - cellPhasePermDiag.push_back(minPerm); - cellPhasePermDiag.push_back(minPerm); - if (helper.upscaleBothPhases) { - cellPhase2PermDiag.push_back(minPerm); - cellPhase2PermDiag.push_back(minPerm); - cellPhase2PermDiag.push_back(minPerm); - } - - if (helper.satnums[cell_idx] > 0) { // handle "no rock" cells with satnum zero - // cout << endl << "Cell no. " << cell_idx << endl; - double PtestvalueCell; - if (includeGravity) { - PtestvalueCell = Ptestvalue - dP[cell_idx]; - } - else { - PtestvalueCell = Ptestvalue; - } - - if (! helper.anisotropic_input) { - - double Jvalue = sqrt(helper.perms[0][cell_idx] * milliDarcyToSqMetre/helper.poros[cell_idx]) * PtestvalueCell; - //cout << "JvalueCell: " << Jvalue << endl; - double waterSaturationCell - = helper.InvJfunctions[int(helper.satnums[cell_idx])-1].evaluate(Jvalue); - waterVolumeLF += waterSaturationCell * helper.cellPoreVolumes[cell_idx]; - - // Compute cell relative permeability. We use a lower cutoff-value as we - // easily divide by zero here. When water saturation is - // zero, we get 'inf', which is circumvented by the cutoff value. - cellPhasePerm = - helper.Krfunctions[0][0][int(helper.satnums[cell_idx])-1].evaluate(waterSaturationCell) * - helper.perms[0][cell_idx]; - if (helper.upscaleBothPhases) { - cellPhase2Perm = - helper.Krfunctions[0][1][int(helper.satnums[cell_idx])-1].evaluate(waterSaturationCell) * - helper.perms[0][cell_idx]; - } - } - else { - double waterSaturationCell = helper.SwPcfunctions[int(helper.satnums[cell_idx])-1].evaluate(PtestvalueCell); - //cout << PtestvalueCell << "\t" << waterSaturationCell << endl; - waterVolumeLF += waterSaturationCell * helper.cellPoreVolumes[cell_idx]; - - cellPhasePermDiag[0] = helper.Krfunctions[0][0][int(helper.satnums[cell_idx])-1].evaluate(waterSaturationCell) * - helper.perms[0][cell_idx]; - cellPhasePermDiag[1] = helper.Krfunctions[1][0][int(helper.satnums[cell_idx])-1].evaluate(waterSaturationCell) * - helper.perms[1][cell_idx]; - cellPhasePermDiag[2] = helper.Krfunctions[2][0][int(helper.satnums[cell_idx])-1].evaluate(waterSaturationCell) * - helper.perms[2][cell_idx]; - if (helper.upscaleBothPhases) { - cellPhase2PermDiag[0] = helper.Krfunctions[0][1][int(helper.satnums[cell_idx])-1].evaluate(waterSaturationCell) * - helper.perms[0][cell_idx]; - cellPhase2PermDiag[1] = helper.Krfunctions[1][1][int(helper.satnums[cell_idx])-1].evaluate(waterSaturationCell) * - helper.perms[1][cell_idx]; - cellPhase2PermDiag[2] = helper.Krfunctions[2][1][int(helper.satnums[cell_idx])-1].evaluate(waterSaturationCell) * - helper.perms[2][cell_idx]; - } - } - - phasePermValues[cell_idx] = cellPhasePerm; - phasePermValuesDiag[cell_idx] = cellPhasePermDiag; - maxPhasePerm = max(maxPhasePerm, cellPhasePerm); - maxPhasePerm = max(maxPhasePerm, *max_element(cellPhasePermDiag.begin(), - cellPhasePermDiag.end())); - if (helper.upscaleBothPhases) { - phase2PermValues[cell_idx] = cellPhase2Perm; - phase2PermValuesDiag[cell_idx] = cellPhase2PermDiag; - maxPhase2Perm = max(maxPhase2Perm, cellPhase2Perm); - maxPhase2Perm = max(maxPhase2Perm, *max_element(cellPhase2PermDiag.begin(), - cellPhase2PermDiag.end())); - } - } - } - // Now we can determine the smallest permitted permeability we can calculate for - - // We have both a fixed bottom limit, as well as a possible higher limit determined - // by a maximum allowable permeability. - double minPhasePerm = max(maxPhasePerm/maxPermContrast, minPerm); - double minPhase2Perm; - if (helper.upscaleBothPhases) minPhase2Perm = max(maxPhase2Perm/maxPermContrast, minPerm); - - // Now remodel the phase permeabilities obeying minPhasePerm - Matrix cellperm(3,3,nullptr); - zero(cellperm); - for (unsigned int i = 0; i < ecl_idx.size(); ++i) { - unsigned int cell_idx = ecl_idx[i]; - zero(cellperm); - if (! helper.anisotropic_input) { - double cellPhasePerm = max(minPhasePerm, phasePermValues[cell_idx]); - accPhasePerm += cellPhasePerm; - double kval = max(minPhasePerm, cellPhasePerm); - cellperm(0,0) = kval; - cellperm(1,1) = kval; - cellperm(2,2) = kval; - } - else { // anisotropic_input - // Truncate values lower than minPhasePerm upwards. - phasePermValuesDiag[cell_idx][0] = max(minPhasePerm, phasePermValuesDiag[cell_idx][0]); - phasePermValuesDiag[cell_idx][1] = max(minPhasePerm, phasePermValuesDiag[cell_idx][1]); - phasePermValuesDiag[cell_idx][2] = max(minPhasePerm, phasePermValuesDiag[cell_idx][2]); - accPhasePerm += phasePermValuesDiag[cell_idx][0]; // not correct anyway - cellperm(0,0) = phasePermValuesDiag[cell_idx][0]; - cellperm(1,1) = phasePermValuesDiag[cell_idx][1]; - cellperm(2,2) = phasePermValuesDiag[cell_idx][2]; - } - helper.upscaler.setPermeability(i, cellperm); - } - - // Output average phase perm, this is just a reality check so that we are not way off. - //cout << ", Arith. mean phase perm = " << accPhasePerm/float(tesselatedCells) << " mD, "; - - // Call single-phase upscaling code - Matrix phasePermTensor = helper.upscaler.upscaleSinglePhase(); - - // Now upscale phase permeability for phase 2 - Matrix phase2PermTensor; - if (helper.upscaleBothPhases) { - zero(cellperm); - for (unsigned int i = 0; i < ecl_idx.size(); ++i) { - unsigned int cell_idx = ecl_idx[i]; - zero(cellperm); - if (! helper.anisotropic_input) { - double cellPhase2Perm = max(minPhase2Perm, phase2PermValues[cell_idx]); - accPhase2Perm += cellPhase2Perm; - double kval = max(minPhase2Perm, cellPhase2Perm); - cellperm(0,0) = kval; - cellperm(1,1) = kval; - cellperm(2,2) = kval; - } - else { // anisotropic_input - // Truncate values lower than minPhasePerm upwards. - phase2PermValuesDiag[cell_idx][0] = max(minPhase2Perm, phase2PermValuesDiag[cell_idx][0]); - phase2PermValuesDiag[cell_idx][1] = max(minPhase2Perm, phase2PermValuesDiag[cell_idx][1]); - phase2PermValuesDiag[cell_idx][2] = max(minPhase2Perm, phase2PermValuesDiag[cell_idx][2]); - accPhase2Perm += phase2PermValuesDiag[cell_idx][0]; // not correct anyway - cellperm(0,0) = phase2PermValuesDiag[cell_idx][0]; - cellperm(1,1) = phase2PermValuesDiag[cell_idx][1]; - cellperm(2,2) = phase2PermValuesDiag[cell_idx][2]; - } - helper.upscaler.setPermeability(i, cellperm); - } - phase2PermTensor = helper.upscaler.upscaleSinglePhase(); - } - - //cout << phasePermTensor << endl; - - - // Here we recalculate the upscaled water saturation, - // although it is already known when we asked for the - // pressure point to compute for. Nonetheless, we - // recalculate here to avoid any minor roundoff-error and - // interpolation error (this means that the saturation - // points are not perfectly uniformly distributed) - helper.WaterSaturation[pointidx] = waterVolumeLF/helper.poreVolume; - - -#ifdef HAVE_MPI - cout << "Rank " << mpi_rank << ": "; -#endif - cout << Ptestvalue << "\t" << helper.WaterSaturation[pointidx]; - // Store and print phase-perm-result - for (int voigtIdx=0; voigtIdx < helper.tensorElementCount; ++voigtIdx) { - helper.PhasePerm[0][pointidx][voigtIdx] = getVoigtValue(phasePermTensor, voigtIdx); - cout << "\t" << getVoigtValue(phasePermTensor, voigtIdx); - if (helper.upscaleBothPhases){ - helper.PhasePerm[1][pointidx][voigtIdx] = getVoigtValue(phase2PermTensor, voigtIdx); - cout << "\t" << getVoigtValue(phase2PermTensor, voigtIdx); - } - } - cout << endl; - } - } - - clock_t finish_upscale_wallclock = clock(); - timeused_upscale_wallclock = (double(finish_upscale_wallclock)-double(start_upscale_wallclock))/CLOCKS_PER_SEC; - - helper.collectResults(); - - // Average time pr. upscaling point: -#ifdef HAVE_MPI - // Sum the upscaling time used by all processes - double timeused_total; - MPI_Reduce(&timeused_upscale_wallclock, &timeused_total, 1, MPI_DOUBLE, - MPI_SUM, 0, MPI_COMM_WORLD); - double avg_upscaling_time_pr_point = timeused_total/(double)helper.points; - -#else - double avg_upscaling_time_pr_point = timeused_upscale_wallclock / (double)helper.points; -#endif - + double timeused_upscaled_wallclock, avg_upscaling_time_pr_point; + std::tie(timeused_upscaled_wallclock, avg_upscaling_time_pr_point) = + helper.upscalePermeability(options, dP, mpi_rank); /* * Step 8c: Make relperm values from phaseperms diff --git a/opm/upscaling/RelPermUtils.cpp b/opm/upscaling/RelPermUtils.cpp index 04c44da..8cff3bd 100644 --- a/opm/upscaling/RelPermUtils.cpp +++ b/opm/upscaling/RelPermUtils.cpp @@ -713,4 +713,257 @@ void RelPermUpscaleHelper::upscaleCapillaryPressure(std::map + RelPermUpscaleHelper::upscalePermeability(std::map& options, + const std::vector& dP, + int mpi_rank) +{ + const double minPerm = atof(options["minPerm"].c_str()); + const double maxPermContrast = atof(options["maxPermContrast"].c_str()); + const double milliDarcyToSqMetre = + Opm::unit::convert::to(1.0*Opm::prefix::milli*Opm::unit::darcy, + Opm::unit::square(Opm::unit::meter)); + + // Put correct number of zeros in, just to be able to access RelPerm[index] later + WaterSaturation.resize(points); + for (size_t i = 0; i < (upscaleBothPhases?2:1); ++i) + PhasePerm[i].resize(points, std::vector(tensorElementCount)); + + // Make vector of capillary pressure points corresponding to uniformly distribued + // saturation points between Swor and Swir. + MonotCubicInterpolator CapPressureVsWaterSaturation(WaterSaturationVsCapPressure.get_fVector(), + WaterSaturationVsCapPressure.get_xVector()); + + for (int pointidx = 1; pointidx <= points; ++pointidx) { + // pointidx=1 corresponds to Swir, pointidx=points to Swor. + double saturation = Swir + (Swor-Swir)/(points-1)*(pointidx-1); + pressurePoints.push_back(CapPressureVsWaterSaturation.evaluate(saturation)); + } + // Preserve max and min pressures + pressurePoints.front() = Pcmax; + pressurePoints.back() = Pcmin; + + // Fill with zeros initially (in case of non-mpi) + node_vs_pressurepoint.resize(points); + +#if HAVE_MPI + // Distribute work load over mpi nodes. + for (int idx=0; idx < points; ++idx) { + // Ensure master node gets equal or less work than the other nodes, since + // master node also computes single phase perm. + node_vs_pressurepoint[idx] = (mpi_nodecount-1) - idx % mpi_nodecount; + } +#endif + + const std::vector& ecl_idx = upscaler.grid().globalCell(); + clock_t start_upscale_wallclock = clock(); + + double waterVolumeLF; + // Now loop through the vector of capillary pressure points that + // this node should compute. + for (int pointidx = 0; pointidx < points; ++pointidx) { + + // Should "I" (mpi-wise) compute this pressure point? + if (node_vs_pressurepoint[pointidx] == mpi_rank) { + + double Ptestvalue = pressurePoints[pointidx]; + + double accPhasePerm = 0.0; + double accPhase2Perm = 0.0; + + double maxPhasePerm = 0.0; + double maxPhase2Perm = 0.0; + + std::vector phasePermValues, phase2PermValues; + std::vector > phasePermValuesDiag, phase2PermValuesDiag; + phasePermValues.resize(satnums.size()); + phasePermValuesDiag.resize(satnums.size()); + if (upscaleBothPhases) { + phase2PermValues.resize(satnums.size()); + phase2PermValuesDiag.resize(satnums.size()); + } + waterVolumeLF = 0.0; + for (size_t i = 0; i < ecl_idx.size(); ++i) { + unsigned int cell_idx = ecl_idx[i]; + double cellPhasePerm = minPerm; + double cellPhase2Perm = minPerm; + std::vector cellPhasePermDiag, cellPhase2PermDiag; + cellPhasePermDiag.resize(3, minPerm); + if (upscaleBothPhases) + cellPhase2PermDiag.resize(3, minPerm); + + if (satnums[cell_idx] > 0) { // handle "no rock" cells with satnum zero + double PtestvalueCell = Ptestvalue; + if (!dP.empty()) + PtestvalueCell -= dP[cell_idx]; + + if (!anisotropic_input) { + double Jvalue = sqrt(perms[0][cell_idx] * milliDarcyToSqMetre/poros[cell_idx]) * PtestvalueCell; + double WaterSaturationCell + = InvJfunctions[int(satnums[cell_idx])-1].evaluate(Jvalue); + waterVolumeLF += WaterSaturationCell * cellPoreVolumes[cell_idx]; + + // Compute cell relative permeability. We use a lower cutoff-value as we + // easily divide by zero here. When water saturation is + // zero, we get 'inf', which is circumvented by the cutoff value. + cellPhasePerm = + Krfunctions[0][0][int(satnums[cell_idx])-1].evaluate(WaterSaturationCell) * + perms[0][cell_idx]; + if (upscaleBothPhases) { + cellPhase2Perm = + Krfunctions[0][1][int(satnums[cell_idx])-1].evaluate(WaterSaturationCell) * + perms[0][cell_idx]; + } + } + else { + double WaterSaturationCell = SwPcfunctions[int(satnums[cell_idx])-1].evaluate(PtestvalueCell); + waterVolumeLF += WaterSaturationCell * cellPoreVolumes[cell_idx]; + + cellPhasePermDiag[0] = Krfunctions[0][0][int(satnums[cell_idx])-1].evaluate(WaterSaturationCell) * + perms[0][cell_idx]; + cellPhasePermDiag[1] = Krfunctions[1][0][int(satnums[cell_idx])-1].evaluate(WaterSaturationCell) * + perms[1][cell_idx]; + cellPhasePermDiag[2] = Krfunctions[2][0][int(satnums[cell_idx])-1].evaluate(WaterSaturationCell) * + perms[2][cell_idx]; + if (upscaleBothPhases) { + cellPhase2PermDiag[0] = Krfunctions[0][1][int(satnums[cell_idx])-1].evaluate(WaterSaturationCell) * + perms[0][cell_idx]; + cellPhase2PermDiag[1] = Krfunctions[1][1][int(satnums[cell_idx])-1].evaluate(WaterSaturationCell) * + perms[1][cell_idx]; + cellPhase2PermDiag[2] = Krfunctions[2][1][int(satnums[cell_idx])-1].evaluate(WaterSaturationCell) * + perms[2][cell_idx]; + } + } + + phasePermValues[cell_idx] = cellPhasePerm; + phasePermValuesDiag[cell_idx] = cellPhasePermDiag; + maxPhasePerm = std::max(maxPhasePerm, cellPhasePerm); + maxPhasePerm = std::max(maxPhasePerm, *std::max_element(cellPhasePermDiag.begin(), + cellPhasePermDiag.end())); + if (upscaleBothPhases) { + phase2PermValues[cell_idx] = cellPhase2Perm; + phase2PermValuesDiag[cell_idx] = cellPhase2PermDiag; + maxPhase2Perm = std::max(maxPhase2Perm, cellPhase2Perm); + maxPhase2Perm = std::max(maxPhase2Perm, *std::max_element(cellPhase2PermDiag.begin(), + cellPhase2PermDiag.end())); + } + } + } + // Now we can determine the smallest permitted permeability we can calculate for + + // We have both a fixed bottom limit, as well as a possible higher limit determined + // by a maximum allowable permeability. + double minPhasePerm = std::max(maxPhasePerm/maxPermContrast, minPerm); + double minPhase2Perm; + if (upscaleBothPhases) + minPhase2Perm = std::max(maxPhase2Perm/maxPermContrast, minPerm); + + // Now remodel the phase permeabilities obeying minPhasePerm + SinglePhaseUpscaler::permtensor_t cellperm(3,3,nullptr); + zero(cellperm); + for (size_t i = 0; i < ecl_idx.size(); ++i) { + unsigned int cell_idx = ecl_idx[i]; + zero(cellperm); + if (!anisotropic_input) { + double cellPhasePerm = std::max(minPhasePerm, phasePermValues[cell_idx]); + accPhasePerm += cellPhasePerm; + double kval = std::max(minPhasePerm, cellPhasePerm); + cellperm(0,0) = kval; + cellperm(1,1) = kval; + cellperm(2,2) = kval; + } + else { // anisotropic_input + // Truncate values lower than minPhasePerm upwards. + phasePermValuesDiag[cell_idx][0] = std::max(minPhasePerm, phasePermValuesDiag[cell_idx][0]); + phasePermValuesDiag[cell_idx][1] = std::max(minPhasePerm, phasePermValuesDiag[cell_idx][1]); + phasePermValuesDiag[cell_idx][2] = std::max(minPhasePerm, phasePermValuesDiag[cell_idx][2]); + accPhasePerm += phasePermValuesDiag[cell_idx][0]; // not correct anyway + cellperm(0,0) = phasePermValuesDiag[cell_idx][0]; + cellperm(1,1) = phasePermValuesDiag[cell_idx][1]; + cellperm(2,2) = phasePermValuesDiag[cell_idx][2]; + } + upscaler.setPermeability(i, cellperm); + } + + // Output average phase perm, this is just a reality check so that we are not way off. + //cout << ", Arith. mean phase perm = " << accPhasePerm/float(tesselatedCells) << " mD, "; + + // Call single-phase upscaling code + SinglePhaseUpscaler::permtensor_t phasePermTensor = upscaler.upscaleSinglePhase(); + + // Now upscale phase permeability for phase 2 + SinglePhaseUpscaler::permtensor_t phase2PermTensor; + if (upscaleBothPhases) { + zero(cellperm); + for (size_t i = 0; i < ecl_idx.size(); ++i) { + unsigned int cell_idx = ecl_idx[i]; + zero(cellperm); + if (!anisotropic_input) { + double cellPhase2Perm = std::max(minPhase2Perm, phase2PermValues[cell_idx]); + accPhase2Perm += cellPhase2Perm; + double kval = std::max(minPhase2Perm, cellPhase2Perm); + cellperm(0,0) = kval; + cellperm(1,1) = kval; + cellperm(2,2) = kval; + } + else { // anisotropic_input + // Truncate values lower than minPhasePerm upwards. + phase2PermValuesDiag[cell_idx][0] = std::max(minPhase2Perm, phase2PermValuesDiag[cell_idx][0]); + phase2PermValuesDiag[cell_idx][1] = std::max(minPhase2Perm, phase2PermValuesDiag[cell_idx][1]); + phase2PermValuesDiag[cell_idx][2] = std::max(minPhase2Perm, phase2PermValuesDiag[cell_idx][2]); + accPhase2Perm += phase2PermValuesDiag[cell_idx][0]; // not correct anyway + cellperm(0,0) = phase2PermValuesDiag[cell_idx][0]; + cellperm(1,1) = phase2PermValuesDiag[cell_idx][1]; + cellperm(2,2) = phase2PermValuesDiag[cell_idx][2]; + } + upscaler.setPermeability(i, cellperm); + } + phase2PermTensor = upscaler.upscaleSinglePhase(); + } + + // Here we recalculate the upscaled water saturation, + // although it is already known when we asked for the + // pressure point to compute for. Nonetheless, we + // recalculate here to avoid any minor roundoff-error and + // interpolation error (this means that the saturation + // points are not perfectly uniformly distributed) + WaterSaturation[pointidx] = waterVolumeLF/poreVolume; + + +#ifdef HAVE_MPI + std::cout << "Rank " << mpi_rank << ": "; +#endif + std::cout << Ptestvalue << "\t" << WaterSaturation[pointidx]; + // Store and print phase-perm-result + for (int voigtIdx=0; voigtIdx < tensorElementCount; ++voigtIdx) { + PhasePerm[0][pointidx][voigtIdx] = getVoigtValue(phasePermTensor, voigtIdx); + std::cout << "\t" << getVoigtValue(phasePermTensor, voigtIdx); + if (upscaleBothPhases){ + PhasePerm[1][pointidx][voigtIdx] = getVoigtValue(phase2PermTensor, voigtIdx); + std::cout << "\t" << getVoigtValue(phase2PermTensor, voigtIdx); + } + } + std::cout << std::endl; + } + } + + clock_t finish_upscale_wallclock = clock(); + double timeused_upscale_wallclock = (double(finish_upscale_wallclock)-double(start_upscale_wallclock))/CLOCKS_PER_SEC; + + collectResults(); + + // Average time pr. upscaling point: +#ifdef HAVE_MPI + // Sum the upscaling time used by all processes + double timeused_total; + MPI_Reduce(&timeused_upscale_wallclock, &timeused_total, 1, MPI_DOUBLE, + MPI_SUM, 0, MPI_COMM_WORLD); + double avg_upscaling_time_pr_point = timeused_total/(double)points; +#else + double avg_upscaling_time_pr_point = timeused_upscale_wallclock / (double)points; +#endif + + return std::make_tuple(timeused_upscale_wallclock, avg_upscaling_time_pr_point); +} + } diff --git a/opm/upscaling/RelPermUtils.hpp b/opm/upscaling/RelPermUtils.hpp index cb126a8..32afb69 100644 --- a/opm/upscaling/RelPermUtils.hpp +++ b/opm/upscaling/RelPermUtils.hpp @@ -27,6 +27,7 @@ #include #include #include +#include #include namespace Opm { @@ -149,6 +150,16 @@ namespace Opm { //! \details Uses the following options: saturationThreshold void upscaleCapillaryPressure(std::map& options, const std::vector& dP); + + //! \brief Upscale permeabilities. + //! \param[in] options Options structure. + //! \param[in] dP Vector with cell pressure gradients. Pass empty vector for no gravity effects. + //! \param[in] mpi_rank MPI rank of this process. + //! \details Uses the following options: minPerm, maxPermContrast + //! \return Tuple with (total time, time per point). + std::tuple upscalePermeability(std::map& options, + const std::vector& dP, + int mpi_rank); private: //! \brief Perform critical saturation check for a single curve. //! \param[in,out] func Function to check for.