#include "nmCalculationDllPebiSolverTask.h" #include "singlePhaseSolver.h" #include "zxLogInstance.h" #include "nmDataAnalyzeManager.h" #include "nmDataWellBase.h" #include "nmDataVerticalWell.h" #include "nmDataVerticalFracturedWell.h" #include "nmDataHorizontalFracturedWell.h" #include "nmDataReservoir.h" #include "nmDataAttribute.h" #include "nmDataRegion.h" #include "nmDataRegionMark.h" #include "nmDataOutline.h" #include "nmDataFracture.h" #include "nmDataFault.h" #include "nmDataTimeStepSetting.h" #include "nmCalculationPebiGrid.h" #include "nmCalculationUtils.h" #include "nmDataAnalyzeManager.h" #include "nmDataPvtParaForPebi.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include /** * @brief 一口求解井在任务启动时捕获的后处理输入。 * * DLL 模型输入已经保存在场景快照中;这里仅保留生成双对数、半对数和结果井位置 * 仍需要的数据,后台线程不再回到 DataManager 查井对象。 */ struct nmPebiSolverWellInputSnapshot { nmPebiSolverWellInputSnapshot() : m_nFlowSectionIndex(1), m_bRateControlled(false), m_bRealWell(false) { } QString m_sWellCode; ///< 项目内稳定井编码。 QVector m_vecFlowPoints; ///< 包含首个零点的原始产量制度。 QPointF m_oLocation; ///< 捕获时的 Map 坐标。 int m_nFlowSectionIndex; ///< 双对数预处理使用的流量段下标。 bool m_bRateControlled; ///< 是否为有产量的主动生产/注入井。 bool m_bRealWell; ///< false 表示仅占 DLL 槽位的手工裂缝。 }; /** * @brief 一次 PEBI 求解任务独占的完整不可变输入。 */ struct nmPebiSolverInputSnapshot { nmPebiSolverInputSnapshot() : m_nSolverType(nmDataAnalyzeManager::PebiSolverOriginal), m_nOmpThreads(1), m_nIluReuseSteps(1), m_bRequiresGridCalculation(false), m_bGridResultNeedsCommit(false), m_bAutoFitTargetOnly(false) { } nmPebiGridInputSnapshot m_oGridInput; ///< 网格、场景、井顺序和授权路径值快照。 nmPebiGridResult m_oGridResult; ///< 已有网格副本或后台新生成的局部网格。 QVector m_vecWellInputs; ///< 与 DLL 井槽位严格对齐。 QVector m_vecPropertyDataSets; ///< 属性插值值副本。 vtkSmartPointer m_pBaseGrid; ///< 完整成果使用的基础 VTK 网格。 int m_nSolverType; ///< PEBI 求解器实现类型。 int m_nOmpThreads; ///< CPU 加速求解线程数。 int m_nIluReuseSteps; ///< ILU 预条件复用步数。 bool m_bRequiresGridCalculation; ///< 后台是否需先用网格值快照计算一次网格。 bool m_bGridResultNeedsCommit; ///< 完整求解成功后是否需在主线程提交新网格。 bool m_bAutoFitTargetOnly; ///< true 时只构造目标井临时曲线,不提交成果。 }; namespace { bool isFiniteSolverNumber(double value) { #ifdef _MSC_VER return _finite(value) != 0; #else return std::isfinite(value); #endif } // 将启用的数据组插值到全部网格单元中心,并覆盖对应的求解器属性数组. bool applyPropertyInterpolation( HX_NWTM_MODEL_INPUT& modelInput, const QVector& dataSets, const QString& licensePath, QString& errorMessage) { bool hasEnabledDataSet = false; for(int i = 0; i < dataSets.size(); ++i) { if(dataSets[i].useForCalculation) { hasEnabledDataSet = true; break; } } if(!hasEnabledDataSet) { return true; } QVector targetPoints; targetPoints.reserve(static_cast(modelInput.GRID.Trinodexy.size())); for(size_t cellIndex = 0; cellIndex < modelInput.GRID.Trinodexy.size(); ++cellIndex) { const dVec1& cellPosition = modelInput.GRID.Trinodexy[cellIndex]; if(cellPosition.size() < 2) { errorMessage = QString("Grid cell %1 has no valid center coordinate.") .arg(static_cast(cellIndex)); return false; } targetPoints.append(QPointF(cellPosition[0], cellPosition[1])); } bool kApplied = false; bool phiApplied = false; bool hApplied = false; for(int dataSetIndex = 0; dataSetIndex < dataSets.size(); ++dataSetIndex) { const nmPropertyInterpolationDataSet& dataSet = dataSets[dataSetIndex]; if(!dataSet.useForCalculation) { continue; } dVec1* solverValues = NULL; bool* propertyApplied = NULL; if(dataSet.property == "k") { solverValues = &modelInput.Base.k; propertyApplied = &kApplied; } else if(dataSet.property == "phi") { solverValues = &modelInput.Base.phi; propertyApplied = &phiApplied; } else if(dataSet.property == "h") { solverValues = &modelInput.Base.h; propertyApplied = &hApplied; } else { errorMessage = QString("Dataset '%1' has an unknown property.") .arg(dataSet.name); return false; } if(*propertyApplied) { errorMessage = QString( "More than one dataset is enabled for property %1.") .arg(dataSet.property); return false; } QVector measurementPoints; QVector measurementValues; measurementPoints.reserve(dataSet.points.size()); measurementValues.reserve(dataSet.points.size()); for(int pointIndex = 0; pointIndex < dataSet.points.size(); ++pointIndex) { const nmPropertyInterpolationPointData& point = dataSet.points[pointIndex]; measurementPoints.append(QPointF(point.x, point.y)); measurementValues.append(point.value); } QVector interpolationValues; QString calculationError; if(!nmCalculationUtils::calculateKriging( targetPoints, measurementPoints, measurementValues, dataSet.nugget, dataSet.sill, dataSet.range, dataSet.model, licensePath, interpolationValues, &calculationError)) { errorMessage = QString("Dataset '%1': %2") .arg(dataSet.name) .arg(calculationError); return false; } solverValues->resize(interpolationValues.size()); for(int valueIndex = 0; valueIndex < interpolationValues.size(); ++valueIndex) { // 属性插值数据层以mD保存;渗透率在写入PEBI数组时转换为D。 (*solverValues)[valueIndex] = dataSet.property == "k" ? nmCalculationUtils::milliDarcyToDarcy( interpolationValues[valueIndex]) : interpolationValues[valueIndex]; } *propertyApplied = true; } return true; } /** * @brief 用网格输出和捕获的场景值组装一次完整模型输入。 */ bool buildModelInputFromSnapshot( const nmPebiSolverInputSnapshot& oSnapshot, const HX_NWTM_GRID_OUTPUT2& oGridOutput, HX_NWTM_MODEL_INPUT& oModelInput, QString& sErrorMessage) { const nmDataBinaryTools::NM_PEBI_SCENE& oScene = oSnapshot.m_oGridInput.m_oScene; oModelInput = HX_NWTM_MODEL_INPUT(oGridOutput); // 第一步:恢复井槽位、产量制度和井筒参数。场景与网格快照在同一次捕获中 // 生成,因此这些外层数组天然使用同一个 DLL 下标顺序。 oModelInput.T = oScene.solverType; oModelInput.Rate.t = oScene.Rate.t; oModelInput.Rate.qo = oScene.Rate.qo; oModelInput.Rate.qg = oScene.Rate.qg; oModelInput.Rate.qw = oScene.Rate.qw; oModelInput.CS.C = oScene.CS.C; oModelInput.CS.S = oScene.CS.S; // 第二步:恢复完整 PVT 数组。不能在工作线程重新读取 PVT 对象,否则界面 // 切换模型或编辑表格时仍会与求解器发生容器并发访问。 oModelInput.PVT.p = oScene.PVT.p; oModelInput.PVT.pb = oScene.PVT.pb; oModelInput.PVT.Rso = oScene.PVT.Rso; oModelInput.PVT.Bo = oScene.PVT.Bo; oModelInput.PVT.Co = oScene.PVT.Co; oModelInput.PVT.miuo = oScene.PVT.miuo; oModelInput.PVT.rouo = oScene.PVT.rouo; oModelInput.PVT.Rv = oScene.PVT.Rv; oModelInput.PVT.Bg = oScene.PVT.Bg; oModelInput.PVT.Cg = oScene.PVT.Cg; oModelInput.PVT.miug = oScene.PVT.miug; oModelInput.PVT.roug = oScene.PVT.roug; oModelInput.PVT.Z = oScene.PVT.Z; oModelInput.PVT.Rsw = oScene.PVT.Rsw; oModelInput.PVT.Bw = oScene.PVT.Bw; oModelInput.PVT.Cw = oScene.PVT.Cw; oModelInput.PVT.miuw = oScene.PVT.miuw; oModelInput.PVT.rouw = oScene.PVT.rouw; oModelInput.PVT.V = oScene.PVT.V; oModelInput.PVT.k_kinitial = oScene.PVT.k_kinitial; oModelInput.PVT.Cf_Cfinitial = oScene.PVT.Cf_Cfinitial; oModelInput.PVT.So = oScene.PVT.So; oModelInput.PVT.Kro = oScene.PVT.Kro; oModelInput.PVT.Sg = oScene.PVT.Sg; oModelInput.PVT.Krg = oScene.PVT.Krg; oModelInput.PVT.Sw = oScene.PVT.Sw; oModelInput.PVT.Krw = oScene.PVT.Krw; // 第三步:先用储层参考值填满全部网格单元,再按捕获的属性数据组覆盖。 oModelInput.Base.Pi = oScene.Base.Pi; oModelInput.Base.Cti = oScene.Base.Cti; oModelInput.Base.Cf = oScene.Base.Cf; oModelInput.Base.Soi = oScene.Base.Soi; oModelInput.Base.Sgi = oScene.Base.Sgi; oModelInput.Base.Swi = oScene.Base.Swi; oModelInput.Base.d = oScene.Base.d; oModelInput.Base.dt_Min = oScene.Base.dt_Min; oModelInput.Base.dt_Max = oScene.Base.dt_Max; const size_t nCellCount = oGridOutput.Trinodexy.size(); oModelInput.Base.k = dVec1(nCellCount, oScene.Base.k_ref); oModelInput.Base.phi = dVec1(nCellCount, oScene.Base.phi_ref); oModelInput.Base.h = dVec1(nCellCount, oScene.Base.h_ref); return applyPropertyInterpolation( oModelInput, oSnapshot.m_vecPropertyDataSets, oSnapshot.m_oGridInput.m_sLicensePath, sErrorMessage); } bool isReasonableLogLogValue(double value) { const double maxReasonableAbsValue = 1.0e12; return isFiniteSolverNumber(value) && value >= -maxReasonableAbsValue && value <= maxReasonableAbsValue; } bool isValidLogLogPoint(const Point& pt) { return isReasonableLogLogValue(pt.x) && isReasonableLogLogValue(pt.y) && isReasonableLogLogValue(pt.z) && pt.x > 0.0 && pt.z >= DBL_EPSILON; } bool isValidSemiLogPoint(const Point& pt) { return isReasonableLogLogValue(pt.x) && pt.x > 0.0 && !pt.pointData.empty() && isReasonableLogLogValue(pt.pointData[0]); } typedef bool (*CalPseudoPressure)(double, double&, int); // mAlgPseudo.dll 是主界面拟压力算法所在模块. // 此处直接解析 iAlgPseuCaller::calPS 的导出符号, 复用主界面的同一套转换实现. // 函数参数依次为原始压力、输出拟压力和拟压力分区编号. CalPseudoPressure getPseudoPressureConverter() { // 静态保存模块和函数地址, 避免每个粒子、每个压力点重复加载和解析 DLL. // configPsAbouts() 初始化的数据及模式也保存在同一个 mAlgPseudo.dll 模块中. static HMODULE module = LoadLibraryW(L"mAlgPseudo.dll"); static CalPseudoPressure converter = module ? reinterpret_cast(GetProcAddress( module, "?calPS@iAlgPseuCaller@@SA_NNAANH@Z")) : nullptr; return converter; } } nmCalculationDllPebiSolverTask::nmCalculationDllPebiSolverTask( QString sPostprocessingDir, nmDataAnalyzeManager* pDataManager, const QString& sAutoFitTargetWellName, QObject *parent): QThread(parent), m_sPostprocessingDir(sPostprocessingDir), m_pDataManager(pDataManager != nullptr ? pDataManager : nmDataAnalyzeManager::getCurrentInstance()), m_pInputSnapshot(new nmPebiSolverInputSnapshot()), m_bManagerUseActive(false), m_bInputSnapshotValid(false), m_nGridInputRevision(0), m_nResultInputRevision(0), m_lastRunSucceeded(false), m_nSolveTimeMs(-1), m_nPebiCount(-1), m_dPendingScalarMin(0.0), m_dPendingScalarMax(0.0), m_bPendingFullResultReady(false) { // 第一步:构造后到 run() 结束前阻止所属成果提前释放 DataManager。 if(m_pDataManager != nullptr) { m_pDataManager->beginBackgroundUse(); m_bManagerUseActive = true; } // 第二步:在创建线程内一次性冻结全部值输入。run() 启动后不得再读取 // DataManager、井、储层、PVT 或属性插值容器。 m_bInputSnapshotValid = captureInputSnapshot(sAutoFitTargetWellName); } nmCalculationDllPebiSolverTask::~nmCalculationDllPebiSolverTask() { // Qt 4.8 没有可靠的外部 DLL 协作取消接口。强制 terminate() 可能让 DLL、 // VTK 或 STL 对象停在持锁/半析构状态,因此析构只能等待 run() 正常退出。 if(isRunning()) { wait(); } releaseDataManagerUse(); delete m_pInputSnapshot; m_pInputSnapshot = nullptr; } void nmCalculationDllPebiSolverTask::run() { m_nSolveTimeMs = -1; m_nPebiCount = -1; m_lastRunSucceeded = this->execute(); // 完成信号可能触发成果关闭;必须在发信号前结束 DataManager 使用期。 releaseDataManagerUse(); emit sig_calculateDone(m_lastRunSucceeded); } void nmCalculationDllPebiSolverTask::releaseDataManagerUse() { if(m_bManagerUseActive && m_pDataManager != nullptr) { m_pDataManager->endBackgroundUse(); m_bManagerUseActive = false; } } bool nmCalculationDllPebiSolverTask::captureInputSnapshot( const QString& sAutoFitTargetWellName) { if(m_pDataManager == nullptr || m_pInputSnapshot == nullptr || QThread::currentThread() != m_pDataManager->thread()) { qWarning() << "PEBI solver input must be captured on the DataManager thread."; return false; } nmDataNumericalAnalysisCase* pAnalysisCase = m_pDataManager->getNumericalAnalysisCase(); if(pAnalysisCase == nullptr) { return false; } // 第一步:冻结几何和求解输入版本。提交阶段会再次核对,期间发生任何编辑 // 都只会让旧任务结果作废,不会影响后台正在读取的值快照。 m_nGridInputRevision = pAnalysisCase->getGridInputRevision(); m_nResultInputRevision = pAnalysisCase->getResultInputRevision(); const nmPebiSolverInputSnapshot oEmptySnapshot; *m_pInputSnapshot = oEmptySnapshot; m_pInputSnapshot->m_bAutoFitTargetOnly = !sAutoFitTargetWellName.isEmpty(); if(m_pInputSnapshot->m_bAutoFitTargetOnly) { nmDataWellBase* pTargetWell = m_pDataManager->findWellByName(sAutoFitTargetWellName); if(pTargetWell == nullptr || pTargetWell->getWellCode().isEmpty()) { qWarning() << "Auto-fit target well is unavailable:" << sAutoFitTargetWellName; return false; } m_sAutoFitTargetWellCode = pTargetWell->getWellCode(); } else { m_sAutoFitTargetWellCode.clear(); } // 第二步:捕获网格输入和求解场景。该调用只复制值,不执行 DLL,也不导出 // CSV,因此任务创建不会等待另一个网格 DLL 或进行无关磁盘写入。 nmCalculationPebiGrid* pGridService = nmCalculationPebiGrid::getInstance(); if(pGridService == nullptr || !pGridService->captureInputSnapshot( m_pDataManager, m_pInputSnapshot->m_oGridInput) || m_pInputSnapshot->m_oGridInput.m_nGridInputRevision != m_nGridInputRevision) { return false; } const QVector& vecWellOrder = m_pInputSnapshot->m_oGridInput.m_vecSolverWellOrder; m_pInputSnapshot->m_vecWellInputs.clear(); m_pInputSnapshot->m_vecWellInputs.reserve(vecWellOrder.size()); // 第三步:捕获结果后处理仍需使用的井数据。数组与网格快照中的求解器顺序 // 一一对应,手工裂缝保留空占位,真实井必须能通过 WellCode 唯一解析。 for(int nWellIndex = 0; nWellIndex < vecWellOrder.size(); ++nWellIndex) { const nmSolverWellRef& oWellRef = vecWellOrder[nWellIndex]; nmPebiSolverWellInputSnapshot oWellInput; oWellInput.m_sWellCode = oWellRef.m_sWellCode; if(oWellRef.m_eEntryKind == NM_SolverEntry_ManualFracture) { m_pInputSnapshot->m_vecWellInputs.append(oWellInput); continue; } nmDataWellBase* pWellData = m_pDataManager->findWellByCode(oWellRef.m_sWellCode); if(pWellData == nullptr) { qWarning() << "Solver well code is missing from Map:" << oWellRef.m_sWellCode; return false; } oWellInput.m_bRealWell = true; oWellInput.m_vecFlowPoints = pWellData->getFlowPoints(); oWellInput.m_nFlowSectionIndex = pWellData->getIndexF(); oWellInput.m_oLocation = QPointF( pWellData->getX().getValue().toDouble(), pWellData->getY().getValue().toDouble()); oWellInput.m_bRateControlled = m_pDataManager->getCalculationWellMode( oWellRef.m_sWellCode) == NM_CaseWell_RateControlled; if(oWellInput.m_bRateControlled && oWellInput.m_vecFlowPoints.size() < 2) { qWarning() << "Rate-controlled well has no valid rate schedule:" << oWellRef.m_sWellCode; return false; } m_pInputSnapshot->m_vecWellInputs.append(oWellInput); } // 第四步:复制属性插值及 DLL 求解配置。这些设置只影响模型求解结果, // 不应在工作线程中再次从 DataManager 查询。 m_pInputSnapshot->m_vecPropertyDataSets = m_pDataManager->getPropertyInterpolationDataSets(); m_pInputSnapshot->m_nSolverType = m_pDataManager->getPebiSolverType(); m_pInputSnapshot->m_nOmpThreads = m_pDataManager->getPebiOmpThreads(); m_pInputSnapshot->m_nIluReuseSteps = m_pDataManager->getPebiIluReuseSteps(); // 第五步:已有有效网格直接复制 DLL 输出;否则只登记“需要计算”,真正的 // 网格 DLL 在后台使用上面的值快照执行一次,不再读写 DataManager。 HX_NWTM_GRID_OUTPUT1 oGridOutput1; HX_NWTM_GRID_OUTPUT2 oGridOutput2; int nPebiCount = -1; if(pGridService->copyCurrentGridFor(m_pDataManager, oGridOutput1, oGridOutput2, nPebiCount)) { m_pInputSnapshot->m_oGridResult.m_oGridOutput1 = oGridOutput1; m_pInputSnapshot->m_oGridResult.m_oGridOutput2 = oGridOutput2; m_pInputSnapshot->m_oGridResult.m_nPebiCount = nPebiCount; m_pInputSnapshot->m_oGridResult.m_bSucceeded = true; if(!m_pInputSnapshot->m_bAutoFitTargetOnly) { m_pInputSnapshot->m_pBaseGrid = m_pDataManager->getUnstructuredGridCopy(); if(m_pInputSnapshot->m_pBaseGrid == nullptr || m_pInputSnapshot->m_pBaseGrid->GetNumberOfCells() <= 0) { return false; } } } else { m_pInputSnapshot->m_bRequiresGridCalculation = true; } return m_pInputSnapshot->m_vecWellInputs.size() == vecWellOrder.size(); } bool nmCalculationDllPebiSolverTask::wasSuccessful() const { return m_lastRunSucceeded; } int nmCalculationDllPebiSolverTask::getSolveTimeMs() const { return m_nSolveTimeMs; } int nmCalculationDllPebiSolverTask::getPebiCount() const { return m_nPebiCount; } QVector > nmCalculationDllPebiSolverTask::getAutoFitResultPressure() const { return m_autoFitResultPressure; } QVector > nmCalculationDllPebiSolverTask::getAutoFitResultLogLog() const { return m_autoFitResultLogLog; } QVector > nmCalculationDllPebiSolverTask::getAutoFitResultSemiLog() const { return m_autoFitResultSemiLog; } bool nmCalculationDllPebiSolverTask::execute() { return this->execPebiMode(); } bool nmCalculationDllPebiSolverTask::execPebiMode() { if(!m_bInputSnapshotValid || m_pInputSnapshot == nullptr) { return false; } // 第一步:清空本次任务自己的输出。旧成果仍保留在 DataManager 中,只有主线程 // 最终提交成功后才会被整体替换。 m_bPendingFullResultReady = false; m_vecPendingWellResults.clear(); m_mapPendingTimeSteps.clear(); m_pPendingResultGrid = nullptr; m_autoFitResultPressure.clear(); m_autoFitResultLogLog.clear(); m_autoFitResultSemiLog.clear(); nmPebiGridResult& oGridResult = m_pInputSnapshot->m_oGridResult; if(m_pInputSnapshot->m_bRequiresGridCalculation) { // 第二步:网格无效时只使用构造阶段捕获的值快照调用一次网格 DLL。 // 自动拟合候选不创建 VTK;完整求解创建局部 VTK,完成后回主线程提交。 const bool bCreateUnstructuredGrid = !m_pInputSnapshot->m_bAutoFitTargetOnly; if(!nmCalculationPebiGrid::getInstance()->calculateSnapshot( m_pInputSnapshot->m_oGridInput, oGridResult, bCreateUnstructuredGrid)) { return false; } if(bCreateUnstructuredGrid) { if(oGridResult.m_pUnstructuredGrid == nullptr || oGridResult.m_pUnstructuredGrid->GetNumberOfCells() <= 0) { return false; } m_pInputSnapshot->m_pBaseGrid = vtkSmartPointer::New(); m_pInputSnapshot->m_pBaseGrid->DeepCopy( oGridResult.m_pUnstructuredGrid); m_pInputSnapshot->m_bGridResultNeedsCommit = true; } } if(!oGridResult.m_bSucceeded || oGridResult.m_oGridOutput1.PEBI_cell.p.empty()) { return false; } // 第三步:由场景值快照重建完整模型输入。属性插值可以耗时,但它只读取任务 // 自己的点集副本,因此放在工作线程不会阻塞或竞争界面数据。 HX_NWTM_MODEL_INPUT oModelInput; QString sInterpolationError; if(!buildModelInputFromSnapshot(*m_pInputSnapshot, oGridResult.m_oGridOutput2, oModelInput, sInterpolationError)) { const QString sLogMessage = QString("Property interpolation failed: %1") .arg(sInterpolationError); qWarning() << sLogMessage; zxLogInstance::getInstance()->writeLogF(sLogMessage); return false; } // 第四步:建网、模型求解和 Kriging 共用 DLL 全局状态,整个 DLL 调用必须串行。 QMutexLocker oDllLocker( nmCalculationUtils::getHxNwtmDllMutex()); HMODULE hModelModule = LoadLibrary(L"HX_NWTM.dll"); if(hModelModule == nullptr) { qWarning() << "Failed to load HX_NWTM.dll. Error code:" << GetLastError(); return false; } typedef void (*HX_NWTM_MODEL_Fun)( HX_NWTM_MODEL_OUTPUT&, const HX_NWTM_MODEL_INPUT&, std::string); typedef void (*SetIntValueFunc)(int); typedef int (*GetIntValueFunc)(); HX_NWTM_MODEL_Fun pfnModel = reinterpret_cast( GetProcAddress(hModelModule, "HX_NWTM_MODEL")); SetIntValueFunc pfnSetSolverType = reinterpret_cast( GetProcAddress(hModelModule, "set_solvetype")); SetIntValueFunc pfnSetOmpThreads = reinterpret_cast( GetProcAddress(hModelModule, "set_omp_threads")); SetIntValueFunc pfnSetIluReuseSteps = reinterpret_cast( GetProcAddress(hModelModule, "set_ilu_reuse_steps")); GetIntValueFunc pfnGetSolveTime = reinterpret_cast( GetProcAddress(hModelModule, "getsolvetime")); if(pfnModel == nullptr || pfnSetSolverType == nullptr || pfnSetOmpThreads == nullptr || pfnSetIluReuseSteps == nullptr || pfnGetSolveTime == nullptr) { qWarning() << "Failed to resolve PEBI solver configuration interface."; FreeLibrary(hModelModule); return false; } HX_NWTM_MODEL_OUTPUT oModelOutput; try { pfnSetSolverType(m_pInputSnapshot->m_nSolverType); if(m_pInputSnapshot->m_nSolverType == nmDataAnalyzeManager::PebiSolverCpuAccelerated) { pfnSetOmpThreads(m_pInputSnapshot->m_nOmpThreads); pfnSetIluReuseSteps(m_pInputSnapshot->m_nIluReuseSteps); } pfnModel(oModelOutput, oModelInput, m_pInputSnapshot->m_oGridInput.m_sLicensePath .toStdString()); m_nPebiCount = oGridResult.m_nPebiCount; m_nSolveTimeMs = pfnGetSolveTime(); } catch(const std::exception& e) { qWarning() << QString("C++ Exception during HX_NWTM_MODEL call: %1") .arg(e.what()); logHX_NWTM_MODEL_INPUT_Simplified(oModelInput); FreeLibrary(hModelModule); return false; } catch(...) { qWarning() << "SEH Exception Occurred during HX_NWTM_MODEL call"; logHX_NWTM_MODEL_INPUT_Simplified(oModelInput); FreeLibrary(hModelModule); return false; } // 后处理只读取本次局部输出,无需继续占用进程级 DLL 锁。 FreeLibrary(hModelModule); oDllLocker.unlock(); // 第五步:把井曲线和场压力构造成任务局部结果。该函数同样只读取输入快照。 const bool bSucceeded = buildPebiModeResult( oModelOutput, oModelInput.T, oGridResult.m_oGridOutput1); return bSucceeded; } std::vector HX_logderivative(const std::vector& x, const std::vector& y, const int& n) { // 功能: 对数导数函数 // 作者: 何辉 // 日期: 2024.07.16 // 单位: 西安华线石油科技有限公司(西安石油大学) std::vector d; d.resize(n - 1); d[0] = 0.5 * (y[1] - y[0]) / (x[1] - x[0]) * (x[1] + x[0]); d[0] = max(d[0], DBL_EPSILON); for (int i = 1; i < n - 2; ++i) { d[i] = (y[i - 1] * (x[i] - x[i + 1]) / ((x[i - 1] - x[i]) * (x[i - 1] - x[i + 1])) + y[i] * (2 * x[i] - x[i - 1] - x[i + 1]) / ((x[i] - x[i - 1]) * (x[i] - x[i + 1])) + y[i + 1] * (x[i] - x[i - 1]) / ((x[i - 1] - x[i + 1]) * (x[i] - x[i + 1]))) * x[i]; d[i] = max(d[i], DBL_EPSILON); } d[n - 2] = 0.5 * (y[n - 2] - y[n - 3]) / (x[n - 2] - x[n - 3]) * (x[n - 2] + x[n - 3]); d[n - 2] = max(d[n - 2], DBL_EPSILON); return d; } bool nmCalculationDllPebiSolverTask::buildPebiModeResult( HX_NWTM_MODEL_OUTPUT& p1, int modelType, const HX_NWTM_GRID_OUTPUT1& oGridOutput) { if(m_pInputSnapshot == nullptr) { return false; } const bool autoFitTargetOnly = m_pInputSnapshot->m_bAutoFitTargetOnly; bool autoFitTargetFound = false; m_autoFitResultPressure.clear(); m_autoFitResultLogLog.clear(); m_autoFitResultSemiLog.clear(); QVector> vvecPressure; QVector> vvecLogLog; QVector> vvecSemiLog; // p1.pw 是求解器返回的原始井底压力, 后续仍按 MPa 保存到 vvecPressure. // 只有气单相变化 PVT 模型的双对数和半对数计算需要改用拟压力. const bool usePseudoPressure = (modelType == static_cast(SMT_Gas_VariablePvt)); CalPseudoPressure calPseudoPressure = usePseudoPressure ? getPseudoPressureConverter() : nullptr; if(usePseudoPressure && calPseudoPressure == nullptr) { qWarning() << "Failed to load the gas pseudo-pressure converter."; return false; } // 井顺序与后处理输入都来自任务启动时的同一份快照。 const QVector& vecWellsOrder = m_pInputSnapshot->m_oGridInput.m_vecSolverWellOrder; const QVector& vecWellInputs = m_pInputSnapshot->m_vecWellInputs; // 第一步:真实井必须逐口具有与公共时间轴等长的井底压力。 // 观察井虽然没有源汇项,也必须由 DLL 返回压力,否则不能形成有效多井结果。 if(vecWellsOrder.isEmpty() || vecWellInputs.size() != vecWellsOrder.size() || p1.t.empty()) { return false; } for(size_t nTimeIndex = 0; nTimeIndex < p1.t.size(); ++nTimeIndex) { if(!isFiniteSolverNumber(p1.t[nTimeIndex])) { qWarning() << "PEBI returned a non-finite time at index:" << static_cast(nTimeIndex); return false; } } for(int nIndex = 0; nIndex < vecWellsOrder.size(); ++nIndex) { const nmSolverWellRef& oWellRef = vecWellsOrder[nIndex]; if(oWellRef.m_eEntryKind == NM_SolverEntry_ManualFracture) { continue; } if(nIndex >= static_cast(p1.pw.size()) || p1.pw[nIndex].empty() || p1.pw[nIndex].size() != p1.t.size() || !vecWellInputs[nIndex].m_bRealWell || vecWellInputs[nIndex].m_sWellCode != oWellRef.m_sWellCode) { qWarning() << "Incomplete well pressure returned for WellCode:" << oWellRef.m_sWellCode; return false; } for(size_t nTimeIndex = 0; nTimeIndex < p1.pw[nIndex].size(); ++nTimeIndex) { if(!isFiniteSolverNumber(p1.pw[nIndex][nTimeIndex])) { qWarning() << "PEBI returned a non-finite pressure for WellCode:" << oWellRef.m_sWellCode << "time index:" << static_cast(nTimeIndex); return false; } } } // 遍历每口井,处理其数据 for(int wellIdx = 0; wellIdx < vecWellsOrder.size(); ++wellIdx) { const nmSolverWellRef& oWellRef = vecWellsOrder[wellIdx]; const nmPebiSolverWellInputSnapshot& oWellInput = vecWellInputs[wellIdx]; if(autoFitTargetOnly && oWellRef.m_sWellCode != m_sAutoFitTargetWellCode) { continue; } // 手工裂缝不是结果井,只用于保持 DLL 下标与网格输入一致。 if(oWellRef.m_eEntryKind == NM_SolverEntry_ManualFracture) { continue; } // 3.1 填充井底压力数据到局部变量 QVector currentWellTime; QVector currentWellPressure; for(size_t i = 0; i < p1.pw[wellIdx].size(); ++i) { currentWellTime.append(p1.t[i]); currentWellPressure.append(p1.pw[wellIdx][i]); } vvecPressure.clear(); // 清空上次循环的数据,为当前井准备 vvecPressure.append(currentWellTime); vvecPressure.append(currentWellPressure); // 3.2 计算双对数和半对数曲线数据。井对象可能已被界面修改,后台只允许 // 使用构造阶段复制出的流量、流动段和坐标。 if(!oWellInput.m_bRealWell || oWellInput.m_sWellCode != oWellRef.m_sWellCode) { return false; } const bool bRateControlled = oWellInput.m_bRateControlled && oWellInput.m_vecFlowPoints.size() >= 2; vvecLogLog.clear(); vvecSemiLog.clear(); if(bRateControlled) { // 准备压力数据 (用于传递给外部 DLL) std::vector wellPressureDataForDll; for(size_t i = 0; i < p1.pw[wellIdx].size(); ++i) { Point pt; pt.x = p1.t[i]; // 使用局部副本转换, 不修改 p1.pw 中需要保存和显示的原始压力. double pressureForLog = p1.pw[wellIdx][i]; if(usePseudoPressure && (!isFiniteSolverNumber(pressureForLog) // -1 与主界面调用一致, 表示使用 configPsAbouts() 配置的当前模式. || !calPseudoPressure(pressureForLog, pressureForLog, -1) || !isFiniteSolverNumber(pressureForLog))) { qWarning() << "Failed to convert gas pressure to pseudo-pressure:" << p1.pw[wellIdx][i] << "well code:" << oWellRef.m_sWellCode; return false; } pt.y = pressureForLog; pt.z = 0.0; wellPressureDataForDll.push_back(pt); } // 准备流量段数据 const QVector& vecTimeQ = oWellInput.m_vecFlowPoints; // 防御空流量数据,避免 -1 转成 std::vector 的巨大无符号长度。 int nTimeNumQ = qMax(0, vecTimeQ.size() - 1); // 移除第一个 0 点 std::vector timeQ(nTimeNumQ); std::vector q(nTimeNumQ); for(int i = 0; i < nTimeNumQ; ++i) { timeQ[i] = vecTimeQ[i + 1].x(); q[i] = vecTimeQ[i + 1].y(); } // 调用外部 DLL 计算双对数曲线 std::vector logPreResultFromDll; // 存储 DLL 的计算结果 const int iSectionFlowIndex = oWellInput.m_nFlowSectionIndex; // 第一步:无产量观察井只接收压力结果,不调用依赖产量制度的曲线 DLL。 HMODULE hMod_solver = nTimeNumQ > 0 ? LoadLibrary(L"singlePhaseSolverDll.dll") : nullptr; if(nTimeNumQ <= 0) { // 第二步:保持固定的数据列结构,便于结果保存和后续读取; // 各列为空明确表示该观察井没有可展示的双对数、半对数结果。 vvecLogLog.clear(); vvecLogLog.append(QVector()); vvecLogLog.append(QVector()); vvecLogLog.append(QVector()); vvecSemiLog.clear(); vvecSemiLog.append(QVector()); vvecSemiLog.append(QVector()); } else if(hMod_solver) { typedef bool (*PreLog)(const std::vector&, const int&, double*, double*, int, std::vector&); PreLog preLogFun = (PreLog)GetProcAddress(hMod_solver, "logLogPre"); if(nullptr == preLogFun) { FreeLibrary(hMod_solver); std::cout << "preLogFun failed!\n"; return false; } // 气井传入拟压力序列, 油井和水井仍传入原始压力序列. // 计算失败必须向上返回, 避免 PSO 使用空曲线或上一粒子的旧曲线. if(!preLogFun(wellPressureDataForDll, iSectionFlowIndex, timeQ.data(), q.data(), nTimeNumQ, logPreResultFromDll)) { FreeLibrary(hMod_solver); return false; } // 填充双对数曲线数据到局部变量 QVector logX, logY, logZ; // 检查结果是否为空,并跳过第一个点 if (logPreResultFromDll.size() > 1) { // 从索引 1 开始遍历,跳过索引 0 的第一个点 for (size_t i = 1; i < logPreResultFromDll.size(); ++i) { const auto& pt = logPreResultFromDll[i]; if(!isValidLogLogPoint(pt)) { continue; } logX.append(pt.x); logY.append(pt.y); logZ.append(pt.z); } } vvecLogLog.clear(); // 清空上次循环的数据 vvecLogLog.append(logX); vvecLogLog.append(logY); vvecLogLog.append(logZ); // 填充半对数曲线数据到局部变量 QVector semiLogX, semiLogY; // 检查结果是否为空,并跳过第一个点 if (logPreResultFromDll.size() > 1) { // 半对数曲线也应该同步跳过第一个点 for (size_t i = 1; i < logPreResultFromDll.size(); ++i) { const auto& pt = logPreResultFromDll[i]; if(!isValidSemiLogPoint(pt)) { continue; } semiLogX.append(pt.x); semiLogY.append(pt.pointData[0]); } } vvecSemiLog.clear(); // 清空上次循环的数据 vvecSemiLog.append(semiLogX); vvecSemiLog.append(semiLogY); FreeLibrary(hMod_solver); } else { qWarning() << "Failed to load singlePhaseSolverDll.dll."; return false; } } if(autoFitTargetOnly) { // 粒子评价结果保存在任务对象内,避免反复改写 DataManager. m_autoFitResultPressure = vvecPressure; m_autoFitResultLogLog = vvecLogLog; m_autoFitResultSemiLog = vvecSemiLog; autoFitTargetFound = true; break; } else { // 所有真实井先写入任务局部快照;观察井的双对数和半对数保持为空。 // 任意后续步骤失败时,整个快照会被丢弃,旧成果不会被局部覆盖。 nmPebiWellResultSnapshot oWellResult; oWellResult.m_sWellCode = oWellRef.m_sWellCode; oWellResult.m_vecPressure = vvecPressure; oWellResult.m_vecLogLog = vvecLogLog; oWellResult.m_vecSemiLog = vvecSemiLog; oWellResult.m_oLocation = oWellInput.m_oLocation; m_vecPendingWellResults.append(oWellResult); } } if(autoFitTargetOnly) { const bool pressureValid = autoFitTargetFound && m_autoFitResultPressure.size() >= 2 && !m_autoFitResultPressure[0].isEmpty() && m_autoFitResultPressure[0].size() == m_autoFitResultPressure[1].size(); const bool logLogValid = m_autoFitResultLogLog.size() >= 3 && !m_autoFitResultLogLog[0].isEmpty() && m_autoFitResultLogLog[0].size() == m_autoFitResultLogLog[1].size() && m_autoFitResultLogLog[0].size() == m_autoFitResultLogLog[2].size(); return pressureValid && logLogValid; } // 第四步:完整求解还必须具有与时间轴、网格绘图单元严格对应的场压力。 if(p1.p.size() != p1.t.size() || oGridOutput.PEBI_cell.isplot.empty()) { return false; } // 计算有效单元数量 size_t actualPlotCellsCount = 0; for(size_t i = 0; i < oGridOutput.PEBI_cell.isplot.size(); ++i) { if(oGridOutput.PEBI_cell.isplot[i] == 1) { actualPlotCellsCount++; } } if(actualPlotCellsCount == 0) { return false; } for(size_t nTimeIndex = 0; nTimeIndex < p1.p.size(); ++nTimeIndex) { if(p1.p[nTimeIndex].size() < oGridOutput.PEBI_cell.isplot.size()) { qWarning() << "Incomplete field pressure returned for time index:" << static_cast(nTimeIndex); return false; } for(size_t nCellIndex = 0; nCellIndex < oGridOutput.PEBI_cell.isplot.size(); ++nCellIndex) { if(oGridOutput.PEBI_cell.isplot[nCellIndex] == 1 && !isFiniteSolverNumber(p1.p[nTimeIndex][nCellIndex])) { qWarning() << "PEBI returned a non-finite field pressure at time/cell:" << static_cast(nTimeIndex) << static_cast(nCellIndex); return false; } } } // 第五步:在局部时间步映射中构造场压力,同时计算全局标量范围。 // 此处不触碰 DataManager,较短的新时间轴会在主线程整体替换时自然清除旧帧。 m_mapPendingTimeSteps.clear(); double dMinP = DBL_MAX; double dMaxP = -DBL_MAX; // 为每个时间步生成完整的 VTK 压力数组。 for(size_t timeIdx = 0; timeIdx < p1.p.size(); ++timeIdx) { // 获取当前时间步的时间值 double currentTime = p1.t[timeIdx]; // 创建 vtkDoubleArray 来存储当前时间步的压力数据 vtkSmartPointer pressureData = vtkSmartPointer::New(); pressureData->SetName("p"); // 设置标量数据的名称 // 直接使用外面计算好的 actualPlotCellsCount pressureData->SetNumberOfValues(actualPlotCellsCount); size_t destIdx = 0; // 目标数组的索引 double dCurrentTimeMin = DBL_MAX; // 当前时间步的最小值 double dCurrentTimeMax = -DBL_MAX; // 当前时间步的最大值 // 遍历原始数据,并根据 isplot 填充到 pressureData for(size_t i = 0; i < oGridOutput.PEBI_cell.isplot.size(); ++i) { if(oGridOutput.PEBI_cell.isplot[i] == 1) { // 只有当 isplot 为1时才考虑这个单元格 pressureData->SetValue(destIdx, p1.p[timeIdx][i]); // 更新当前时间步的范围 dCurrentTimeMin = qMin(dCurrentTimeMin, p1.p[timeIdx][i]); dCurrentTimeMax = qMax(dCurrentTimeMax, p1.p[timeIdx][i]); destIdx++; } } // 更新全局范围 dMinP = qMin(dMinP, dCurrentTimeMin); dMaxP = qMax(dMaxP, dCurrentTimeMax); m_mapPendingTimeSteps.insert(currentTime, pressureData); } // QMap 使用时间作为唯一键;重复时间会覆盖旧帧,因此必须检查数量一致。 if(m_mapPendingTimeSteps.size() != static_cast(p1.t.size())) { return false; } // 第六步:基础网格在任务创建或局部建网阶段已经完成深拷贝。后台结果构造 // 不再读取 DataManager 中可能被网格窗口替换的 VTK 指针。 m_pPendingResultGrid = m_pInputSnapshot->m_pBaseGrid; if(m_pPendingResultGrid == nullptr || m_pPendingResultGrid->GetNumberOfCells() <= 0) { return false; } int nExpectedWellResultCount = 0; for(int nIndex = 0; nIndex < vecWellsOrder.size(); ++nIndex) { if(vecWellsOrder[nIndex].m_eEntryKind == NM_SolverEntry_Well) { ++nExpectedWellResultCount; } } if(m_vecPendingWellResults.size() != nExpectedWellResultCount) { return false; } m_dPendingScalarMin = dMinP; m_dPendingScalarMax = dMaxP; m_bPendingFullResultReady = true; return true; } bool nmCalculationDllPebiSolverTask::commitResult( nmDataAnalyzeManager* pDataManager) { // 第一步:自动拟合只返回目标曲线,不允许提交完整成果;手工求解必须回到 // DataManager 所属线程执行,保证界面看不到逐项替换过程中的中间状态。 if(m_pInputSnapshot == nullptr || m_pInputSnapshot->m_bAutoFitTargetOnly || !m_bPendingFullResultReady || pDataManager == nullptr || pDataManager != m_pDataManager || QThread::currentThread() != pDataManager->thread()) { return false; } nmDataNumericalAnalysisCase* pAnalysisCase = pDataManager->getNumericalAnalysisCase(); if(pAnalysisCase == nullptr || pAnalysisCase->getGridInputRevision() != m_nGridInputRevision || pAnalysisCase->getResultInputRevision() != m_nResultInputRevision || m_mapPendingTimeSteps.isEmpty() || m_pPendingResultGrid == nullptr || m_pPendingResultGrid->GetNumberOfCells() <= 0) { return false; } // 第二步:在改动旧成果前一次性解析全部目标井,并检查 WellCode 不重复。 QVector vecTargetWells; QSet setWellCodes; vecTargetWells.reserve(m_vecPendingWellResults.size()); for(int nIndex = 0; nIndex < m_vecPendingWellResults.size(); ++nIndex) { const nmPebiWellResultSnapshot& oWellResult = m_vecPendingWellResults[nIndex]; nmDataWellBase* pWellData = pDataManager->findWellByCode(oWellResult.m_sWellCode); if(pWellData == nullptr || setWellCodes.contains(oWellResult.m_sWellCode)) { return false; } setWellCodes.insert(oWellResult.m_sWellCode); vecTargetWells.append(pWellData); } // 第三步:任务后台生成了新网格时,先在当前主线程按同一输入版本提交。 // 已有网格路径则再次确认其仍然有效。两条路径都不允许旧任务覆盖新编辑。 if(m_pInputSnapshot->m_bGridResultNeedsCommit) { if(!nmCalculationPebiGrid::getInstance()->commitSnapshotResult( pDataManager, m_pInputSnapshot->m_oGridInput, m_pInputSnapshot->m_oGridResult)) { return false; } } else if(!pAnalysisCase->isGridValid()) { return false; } // 第四步:最后一次登记版本。登记失败时旧成果仍未被修改;登记成功后当前 // 主线程事件不会被其他编辑操作插入,因此后续替换不存在可恢复失败分支。 if(!pAnalysisCase->markResultsAvailableIfCurrent( m_nGridInputRevision, m_nResultInputRevision)) { return false; } // 第五步:版本检查和对象解析全部通过后,在当前主线程事件内整体替换。 // 这些 setter 不发事件也不包含可恢复失败分支,外部只能在本函数返回后看到新成果。 pDataManager->clearWellLocations(); for(int nIndex = 0; nIndex < m_vecPendingWellResults.size(); ++nIndex) { const nmPebiWellResultSnapshot& oWellResult = m_vecPendingWellResults[nIndex]; nmDataWellBase* pWellData = vecTargetWells[nIndex]; pWellData->setResultPressure(oWellResult.m_vecPressure); pWellData->setResultLogLog(oWellResult.m_vecLogLog); pWellData->setResultSemiLog(oWellResult.m_vecSemiLog); pDataManager->addWellLocation(oWellResult.m_sWellCode, oWellResult.m_oLocation); } pDataManager->clearTimeSteps(); QMap >::const_iterator oTimeIt = m_mapPendingTimeSteps.constBegin(); for(; oTimeIt != m_mapPendingTimeSteps.constEnd(); ++oTimeIt) { pDataManager->addTimeStep(oTimeIt.key(), oTimeIt.value()); } pDataManager->setScalarRangeP(m_dPendingScalarMin, m_dPendingScalarMax); pDataManager->setResultBaseGrid(m_pPendingResultGrid); // 第六步:清除待提交标志,防止重复完成信号再次覆盖后续结果。 m_bPendingFullResultReady = false; return true; } //bool nmCalculationDllPebiSolverTask::savePebiModeResult(HX_NWTM_MODEL_OUTPUT& p1) //{ // nmDataAnalyzeManager* pDataInstance = nmDataAnalyzeManager::getCurrentInstance(); // // QVector> vvecPressure; // QVector> vvecLogLog; // QVector> vvecSemiLog; // // // 获取参与求解的井的顺序 // QVector> vecWellsOrder = nmDataAnalyzeManager::getCurrentInstance()->getCalculationWells(); // // // 清空井名和二维位置的映射 // pDataInstance->clearWellLocations(); // // // 遍历每口井,处理其数据 // for(int wellIdx = 0; wellIdx < vecWellsOrder.size(); ++wellIdx) { // NM_WELL_MODEL eWellType = vecWellsOrder[wellIdx].first; // 获取井的类型 // QString sWellName = vecWellsOrder[wellIdx].second; // 获取井的名称 // // // 跳过裂缝(或未知井类型) // if(eWellType == NM_WELL_MODEL::Unknow_Well) { // continue; // } // // // 3.1 填充井底压力数据到局部变量 // QVector currentWellTime; // QVector currentWellPressure; // P_wf(t) // // // 确保 p1.pw[wellIdx] 存在且大小与 p1.t 匹配 // if(wellIdx < p1.pw.size()) { // for(size_t i = 0; i < p1.pw[wellIdx].size(); ++i) { // if(i < p1.t.size()) { // 确保时间数据也存在 // currentWellTime.append(p1.t[i]); // currentWellPressure.append(p1.pw[wellIdx][i]); // } // } // } // // vvecPressure.clear(); // 清空上次循环的数据,为当前井准备 // vvecPressure.append(currentWellTime); // vvecPressure.append(currentWellPressure); // // // 3.2 计算双对数和半对数曲线数据并存储到局部变量 // nmDataWellBase* pWellData = nmDataAnalyzeManager::getCurrentInstance()->findWellByName(sWellName); // // if(pWellData) { // // // --- START: 替换为 HX_logderivative 计算双对数数据 --- // // // 1. 获取初始压力 Pi,用于计算压力降 DeltaP // double initialPressure = nmDataAnalyzeManager::getCurrentInstance()->getReservoirData()->getInitialPressure().getValue().toDouble(); // // // 为了演示,如果无法获取 Pi,我们将使用数组中第一个压力点作为初始压力 Pi // //if (!currentWellPressure.isEmpty()) { // // // 通常第一个压力点 Pwf(t=0) 就是初始压力 Pi // // initialPressure = currentWellPressure.first(); // //} // // // 2. 准备 t (x) 和 DeltaP (y) 的 std::vector // std::vector t_vec; // 时间 t (x) // std::vector deltaP_vec; // 压力降 DeltaP = Pi - Pwf(t) (y) // int N = currentWellTime.size(); // // for (int i = 0; i < N; ++i) { // t_vec.push_back(currentWellTime[i]); // // 计算压力降 // deltaP_vec.push_back(initialPressure - currentWellPressure[i]); // } // // // 3. 调用 HX_logderivative 计算对数导数 (P') // std::vector logDerivative_vec; // // if (N > 2) { // 确保有足够的数据点 (N-1 个导数点) // logDerivative_vec = HX_logderivative(t_vec, deltaP_vec, N); // } else { // // 数据点不足,无法计算导数 // logDerivative_vec.clear(); // } // // // 4. 组装双对数曲线数据 (t, DeltaP, P'),严格遵循绘图规则 // QVector logX, logY, logZ; // logX: t_plot, logY: DeltaP_plot, logZ: P' // size_t M = logDerivative_vec.size(); // 导数点数量 M = N - 1 // // if (M >= 1) { // // // --- 4.1 第一个点 (索引 i=0) --- // // 导数时间: t_plot = 0.5 * (x[1] + x[0]) // logX.append(0.5 * (t_vec[1] + t_vec[0])); // // 压力降 (近似): DeltaP[0] 和 DeltaP[1] 的平均值 // logY.append((deltaP_vec[0] + deltaP_vec[1]) / 2.0); // // 导数: d[0] // logZ.append(logDerivative_vec[0]); // // // // --- 4.2 中间点 (索引 i=1 到 N-3) --- // // 在 logDerivative_vec 中,这些点对应索引 i=1 到 M-1 (即 N-3) // for (size_t i = 1; i < (size_t)N - 2; ++i) { // // 导数时间: t_plot = x[i] // logX.append(t_vec[i]); // // 压力降: 匹配原始数据 DeltaP[i] // logY.append(deltaP_vec[i]); // // 导数: d[i] // logZ.append(logDerivative_vec[i]); // } // // // --- 4.3 最后一个点 (索引 i=N-2) --- // if (M >= 2) { // 确保至少有 2 个导数点 // size_t lastIdx = M - 1; // 对应 logDerivative_vec[N-2] // // 导数时间: t_plot = 0.5 * (x[n-2] + x[n-3]) // // 对应 t_vec[N-2] 和 t_vec[N-3] // logX.append(0.5 * (t_vec[N - 2] + t_vec[N - 3])); // // // 压力降 (近似): DeltaP[N-2] 和 DeltaP[N-3] 的平均值 // logY.append((deltaP_vec[N - 2] + deltaP_vec[N - 3]) / 2.0); // // // 导数: d[N-2] // logZ.append(logDerivative_vec[lastIdx]); // } // } // // // 将组装好的数据存储到 vvecLogLog // vvecLogLog.clear(); // 清空上次循环的数据 // vvecLogLog.append(logX); // 时间 (t_plot) // vvecLogLog.append(logY); // 压力降 (DeltaP_plot) // vvecLogLog.append(logZ); // 对数导数 (P') // // // // 5. 组装半对数曲线数据 // // 半对数图 Y 轴是井底压力 Pwf(t),X 轴是时间 t // // vvecSemiLog 存储 [t, Pwf(t)],使用原始数据即可 // vvecSemiLog.clear(); // vvecSemiLog.append(currentWellTime); // X轴: 时间 t // vvecSemiLog.append(currentWellPressure); // Y轴: 井底压力 Pwf(t) // // // --- END: 替换为 HX_logderivative 计算双对数数据 --- // // } // end if(pWellData) // // // 将计算结果保存到对应的井数据里 // // 压力 // pWellData->setResultPressure(vvecPressure); // // 双对数 // pWellData->setResultLogLog(vvecLogLog); // // 半对数 // pWellData->setResultSemiLog(vvecSemiLog); // // // 存储当前井名称和二维位置到映射 // QPointF ptWellCoords(pWellData->getX().getValue().toDouble(), pWellData->getY().getValue().toDouble()); // pDataInstance->addWellLocation(sWellName, ptWellCoords); // } // 主日志函数:打印 HX_NWTM_MODEL_INPUT 的内容 void nmCalculationDllPebiSolverTask::logHX_NWTM_MODEL_INPUT_Simplified(const HX_NWTM_MODEL_INPUT& p0) { // 1. Basic Parameters (T) qDebug() << QString("Parameter T: %1").arg(p0.T); // 2. GRID Data qDebug() << "\n--- GRID Data ---"; qDebug() << QString("GRID.Trinodexy: Capacity %1").arg(p0.GRID.Trinodexy.size()); qDebug() << QString("GRID.Area: Capacity %1").arg(p0.GRID.Area.size()); qDebug() << QString("GRID.D: Capacity %1").arg(p0.GRID.D.size()); // 2.1 GRID.ZhiJingNeiBianJie qDebug() << "\n -- GRID.ZhiJingNeiBianJie --"; qDebug() << QString(" n: %1").arg(p0.GRID.ZhiJingNeiBianJie.n); qDebug() << QString(" XiLinw: Capacity %1").arg(p0.GRID.ZhiJingNeiBianJie.XiLinw.size()); qDebug() << QString(" lw: Capacity %1").arg(p0.GRID.ZhiJingNeiBianJie.lw.size()); qDebug() << QString(" dw: Capacity %1").arg(p0.GRID.ZhiJingNeiBianJie.dw.size()); qDebug() << QString(" rw: Capacity %1").arg(p0.GRID.ZhiJingNeiBianJie.rw.size()); qDebug() << QString(" inwell: Capacity %1").arg(p0.GRID.ZhiJingNeiBianJie.inwell.size()); // 2.2 GRID.LieFengJingNeiBianJie qDebug() << "\n -- GRID.LieFengJingNeiBianJie --"; qDebug() << QString(" n: %1").arg(p0.GRID.LieFengJingNeiBianJie.n); qDebug() << QString(" XiLinf: Capacity %1").arg(p0.GRID.LieFengJingNeiBianJie.XiLinf.size()); qDebug() << QString(" lf: Capacity %1").arg(p0.GRID.LieFengJingNeiBianJie.lf.size()); qDebug() << QString(" df: Capacity %1").arg(p0.GRID.LieFengJingNeiBianJie.df.size()); qDebug() << QString(" xf: Capacity %1").arg(p0.GRID.LieFengJingNeiBianJie.xf.size()); qDebug() << QString(" infra: Capacity %1").arg(p0.GRID.LieFengJingNeiBianJie.infra.size()); // 2.3 GRID.DuoJiYaLieShuiPingJingNeiBianJie qDebug() << "\n -- GRID.DuoJiYaLieShuiPingJingNeiBianJie --"; qDebug() << QString(" n: %1").arg(p0.GRID.DuoJiYaLieShuiPingJingNeiBianJie.n); qDebug() << QString(" XiLinh: Capacity %1").arg(p0.GRID.DuoJiYaLieShuiPingJingNeiBianJie.XiLinh.size()); qDebug() << QString(" lh: Capacity %1").arg(p0.GRID.DuoJiYaLieShuiPingJingNeiBianJie.lh.size()); qDebug() << QString(" dh: Capacity %1").arg(p0.GRID.DuoJiYaLieShuiPingJingNeiBianJie.dh.size()); qDebug() << QString(" dsxf: Capacity %1").arg(p0.GRID.DuoJiYaLieShuiPingJingNeiBianJie.dsxf.size()); qDebug() << QString(" inhor: Capacity %1").arg(p0.GRID.DuoJiYaLieShuiPingJingNeiBianJie.inhor.size()); // 2.4 GRID.WaiBianJie qDebug() << "\n -- GRID.WaiBianJie --"; qDebug() << QString(" n: %1").arg(p0.GRID.WaiBianJie.n); qDebug() << QString(" WaiBianh: Capacity %1").arg(p0.GRID.WaiBianJie.WaiBianh.size()); qDebug() << QString(" WaiBianl: Capacity %1").arg(p0.GRID.WaiBianJie.WaiBianl.size()); qDebug() << QString(" WaiBiand: Capacity %1").arg(p0.GRID.WaiBianJie.WaiBiand.size()); // 2.5 GRID.NeiBuDuanCeng qDebug() << "\n -- GRID.NeiBuDuanCeng --"; qDebug() << QString(" n: %1").arg(p0.GRID.NeiBuDuanCeng.n); qDebug() << QString(" faultb1: Capacity %1").arg(p0.GRID.NeiBuDuanCeng.faultb1.size()); qDebug() << QString(" faultb2: Capacity %1").arg(p0.GRID.NeiBuDuanCeng.faultb2.size()); qDebug() << QString(" faultl1: Capacity %1").arg(p0.GRID.NeiBuDuanCeng.faultl1.size()); qDebug() << QString(" faultd1: Capacity %1").arg(p0.GRID.NeiBuDuanCeng.faultd1.size()); // 2.6 GRID.YuChuLiJuZhen qDebug() << "\n -- GRID.YuChuLiJuZhen --"; qDebug() << QString(" ia: Capacity %1").arg(p0.GRID.YuChuLiJuZhen.ia.size()); qDebug() << QString(" ja: Capacity %1").arg(p0.GRID.YuChuLiJuZhen.ja.size()); qDebug() << QString(" nzeros: Capacity %1").arg(p0.GRID.YuChuLiJuZhen.nzeros.size()); qDebug() << QString(" numk: %1").arg(p0.GRID.YuChuLiJuZhen.numk); // 3. Rate Data qDebug() << "\n--- Rate Data ---"; qDebug() << QString("Rate.t: Contains %1 data sets").arg(p0.Rate.t.size()); qDebug() << QString("Rate.qo: Contains %1 data sets").arg(p0.Rate.qo.size()); qDebug() << QString("Rate.qg: Contains %1 data sets").arg(p0.Rate.qg.size()); qDebug() << QString("Rate.qw: Contains %1 data sets").arg(p0.Rate.qw.size()); // Optional: Print capacity of each inner vector, e.g.: // for (int i = 0; i < p0.Rate.t.size(); ++i) { // zxLogInstance::getInstance()->writeLogF(QString(" Rate.t[%1] Capacity: %2").arg(i).arg(p0.Rate.t[i].size())); // } // 4. Pressure Data qDebug() << "\n--- Pressure Data ---"; qDebug() << QString("Pressure.t: Contains %1 data sets").arg(p0.Pressure.t.size()); qDebug() << QString("Pressure.p: Contains %1 data sets").arg(p0.Pressure.p.size()); // 5. CS Wellbore Storage and Skin Data qDebug() << "\n--- CS Data ---"; qDebug() << QString("CS.C: Capacity %1").arg(p0.CS.C.size()); qDebug() << QString("CS.S: Capacity %1").arg(p0.CS.S.size()); // 6. PVT Fluid Property Data qDebug() << "\n--- PVT Data ---"; qDebug() << QString("PVT.p: Capacity %1").arg(p0.PVT.p.size()); qDebug() << QString("PVT.Rso: Capacity %1").arg(p0.PVT.Rso.size()); qDebug() << QString("PVT.Bo: Capacity %1").arg(p0.PVT.Bo.size()); qDebug() << QString("PVT.Co: Capacity %1").arg(p0.PVT.Co.size()); qDebug() << QString("PVT.miuo: Capacity %1").arg(p0.PVT.miuo.size()); qDebug() << QString("PVT.rouo: Capacity %1").arg(p0.PVT.rouo.size()); qDebug() << QString("PVT.Rv: Capacity %1").arg(p0.PVT.Rv.size()); qDebug() << QString("PVT.Bg: Capacity %1").arg(p0.PVT.Bg.size()); qDebug() << QString("PVT.Cg: Capacity %1").arg(p0.PVT.Cg.size()); qDebug() << QString("PVT.miug: Capacity %1").arg(p0.PVT.miug.size()); qDebug() << QString("PVT.roug: Capacity %1").arg(p0.PVT.roug.size()); qDebug() << QString("PVT.Z: Capacity %1").arg(p0.PVT.Z.size()); qDebug() << QString("PVT.Rsw: Capacity %1").arg(p0.PVT.Rsw.size()); qDebug() << QString("PVT.Bw: Capacity %1").arg(p0.PVT.Bw.size()); qDebug() << QString("PVT.Cw: Capacity %1").arg(p0.PVT.Cw.size()); qDebug() << QString("PVT.miuw: Capacity %1").arg(p0.PVT.miuw.size()); qDebug() << QString("PVT.rouw: Capacity %1").arg(p0.PVT.rouw.size()); qDebug() << QString("PVT.V: Capacity %1").arg(p0.PVT.V.size()); qDebug() << QString("PVT.k_kinitial: Capacity %1").arg(p0.PVT.k_kinitial.size()); qDebug() << QString("PVT.Cf_Cfinitial: Capacity %1").arg(p0.PVT.Cf_Cfinitial.size()); qDebug() << QString("PVT.So: Capacity %1").arg(p0.PVT.So.size()); qDebug() << QString("PVT.Kro: Capacity %1").arg(p0.PVT.Kro.size()); qDebug() << QString("PVT.Sg: Capacity %1").arg(p0.PVT.Sg.size()); qDebug() << QString("PVT.Krg: Capacity %1").arg(p0.PVT.Krg.size()); qDebug() << QString("PVT.Sw: Capacity %1").arg(p0.PVT.Sw.size()); qDebug() << QString("PVT.Krw: Capacity %1").arg(p0.PVT.Krw.size()); // 7. Base Data qDebug() << "\n--- Base Data ---"; qDebug() << QString("Base.Pi: %1").arg(p0.Base.Pi); qDebug() << QString("Base.Cti: %1").arg(p0.Base.Cti); qDebug() << QString("Base.Cf: %1").arg(p0.Base.Cf); qDebug() << QString("Base.Soi: %1").arg(p0.Base.Soi); qDebug() << QString("Base.Sgi: %1").arg(p0.Base.Sgi); qDebug() << QString("Base.Swi: %1").arg(p0.Base.Swi); qDebug() << QString("Base.k: Capacity %1").arg(p0.Base.k.size()); qDebug() << QString("Base.phi: Capacity %1").arg(p0.Base.phi.size()); qDebug() << QString("Base.h: Capacity %1").arg(p0.Base.h.size()); qDebug() << "\n--- End of HX_NWTM_MODEL_INPUT Content ---"; } bool nmCalculationDllPebiSolverTask::saveHX_NWTM_MODEL_INPUT_ToTxt(const HX_NWTM_MODEL_INPUT& p0, const QString& filePath) { QFile file(filePath); if(!file.open(QIODevice::WriteOnly | QIODevice::Text | QIODevice::Truncate)) { // Log an error if the file cannot be opened qDebug() << QString("ERROR: Could not open file for writing: %1").arg(filePath); return false; } QTextStream out(&file); out << "--- Start Printing HX_NWTM_MODEL_INPUT Content ---\n"; // 1. Basic Parameters (T) out << QString("Parameter T: %1\n").arg(p0.T); // 2. GRID Data out << "\n--- GRID Data ---\n"; out << QString("GRID.Trinodexy: Capacity %1\n").arg(p0.GRID.Trinodexy.size()); out << QString("GRID.Area: Capacity %1\n").arg(p0.GRID.Area.size()); out << QString("GRID.D: Capacity %1\n").arg(p0.GRID.D.size()); // 2.1 GRID.ZhiJingNeiBianJie out << "\n -- GRID.ZhiJingNeiBianJie --\n"; out << QString(" n: %1\n").arg(p0.GRID.ZhiJingNeiBianJie.n); out << QString(" XiLinw: Capacity %1\n").arg(p0.GRID.ZhiJingNeiBianJie.XiLinw.size()); out << QString(" lw: Capacity %1\n").arg(p0.GRID.ZhiJingNeiBianJie.lw.size()); out << QString(" dw: Capacity %1\n").arg(p0.GRID.ZhiJingNeiBianJie.dw.size()); out << QString(" rw: Capacity %1\n").arg(p0.GRID.ZhiJingNeiBianJie.rw.size()); out << QString(" inwell: Capacity %1\n").arg(p0.GRID.ZhiJingNeiBianJie.inwell.size()); // 2.2 GRID.LieFengJingNeiBianJie out << "\n -- GRID.LieFengJingNeiBianJie --\n"; out << QString(" n: %1\n").arg(p0.GRID.LieFengJingNeiBianJie.n); out << QString(" XiLinf: Capacity %1\n").arg(p0.GRID.LieFengJingNeiBianJie.XiLinf.size()); out << QString(" lf: Capacity %1\n").arg(p0.GRID.LieFengJingNeiBianJie.lf.size()); out << QString(" df: Capacity %1\n").arg(p0.GRID.LieFengJingNeiBianJie.df.size()); out << QString(" xf: Capacity %1\n").arg(p0.GRID.LieFengJingNeiBianJie.xf.size()); out << QString(" infra: Capacity %1\n").arg(p0.GRID.LieFengJingNeiBianJie.infra.size()); // 2.3 GRID.DuoJiYaLieShuiPingJingNeiBianJie out << "\n -- GRID.DuoJiYaLieShuiPingJingNeiBianJie --\n"; out << QString(" n: %1\n").arg(p0.GRID.DuoJiYaLieShuiPingJingNeiBianJie.n); out << QString(" XiLinh: Capacity %1\n").arg(p0.GRID.DuoJiYaLieShuiPingJingNeiBianJie.XiLinh.size()); out << QString(" lh: Capacity %1\n").arg(p0.GRID.DuoJiYaLieShuiPingJingNeiBianJie.lh.size()); out << QString(" dh: Capacity %1\n").arg(p0.GRID.DuoJiYaLieShuiPingJingNeiBianJie.dh.size()); out << QString(" dsxf: Capacity %1\n").arg(p0.GRID.DuoJiYaLieShuiPingJingNeiBianJie.dsxf.size()); out << QString(" inhor: Capacity %1\n").arg(p0.GRID.DuoJiYaLieShuiPingJingNeiBianJie.inhor.size()); // 2.4 GRID.WaiBianJie out << "\n -- GRID.WaiBianJie --\n"; out << QString(" n: %1\n").arg(p0.GRID.WaiBianJie.n); out << QString(" WaiBianh: Capacity %1\n").arg(p0.GRID.WaiBianJie.WaiBianh.size()); out << QString(" WaiBianl: Capacity %1\n").arg(p0.GRID.WaiBianJie.WaiBianl.size()); out << QString(" WaiBiand: Capacity %1\n").arg(p0.GRID.WaiBianJie.WaiBiand.size()); // 2.5 GRID.NeiBuDuanCeng out << "\n -- GRID.NeiBuDuanCeng --\n"; out << QString(" n: %1\n").arg(p0.GRID.NeiBuDuanCeng.n); out << QString(" faultb1: Capacity %1\n").arg(p0.GRID.NeiBuDuanCeng.faultb1.size()); out << QString(" faultb2: Capacity %1\n").arg(p0.GRID.NeiBuDuanCeng.faultb2.size()); out << QString(" faultl1: Capacity %1\n").arg(p0.GRID.NeiBuDuanCeng.faultl1.size()); out << QString(" faultd1: Capacity %1\n").arg(p0.GRID.NeiBuDuanCeng.faultd1.size()); // 2.6 GRID.YuChuLiJuZhen out << "\n -- GRID.YuChuLiJuZhen --\n"; out << QString(" ia: Capacity %1\n").arg(p0.GRID.YuChuLiJuZhen.ia.size()); out << QString(" ja: Capacity %1\n").arg(p0.GRID.YuChuLiJuZhen.ja.size()); out << QString(" nzeros: Capacity %1\n").arg(p0.GRID.YuChuLiJuZhen.nzeros.size()); out << QString(" numk: %1\n").arg(p0.GRID.YuChuLiJuZhen.numk); // 3. Rate Data out << "\n--- Rate Data ---\n"; out << QString("Rate.t: Contains %1 data sets\n").arg(p0.Rate.t.size()); out << QString("Rate.qo: Contains %1 data sets\n").arg(p0.Rate.qo.size()); out << QString("Rate.qg: Contains %1 data sets\n").arg(p0.Rate.qg.size()); out << QString("Rate.qw: Contains %1 data sets\n").arg(p0.Rate.qw.size()); // Optional: Print capacity of each inner vector, e.g.: // for (int i = 0; i < p0.Rate.t.size(); ++i) { // out << QString(" Rate.t[%1] Capacity: %2\n").arg(i).arg(p0.Rate.t[i].size()); // } // 4. Pressure Data out << "\n--- Pressure Data ---\n"; out << QString("Pressure.t: Contains %1 data sets\n").arg(p0.Pressure.t.size()); out << QString("Pressure.p: Contains %1 data sets\n").arg(p0.Pressure.p.size()); // 5. CS Wellbore Storage and Skin Data out << "\n--- CS Data ---\n"; out << QString("CS.C: Capacity %1\n").arg(p0.CS.C.size()); out << QString("CS.S: Capacity %1\n").arg(p0.CS.S.size()); // 6. PVT Fluid Property Data out << "\n--- PVT Data ---\n"; out << QString("PVT.p: Capacity %1\n").arg(p0.PVT.p.size()); out << QString("PVT.Rso: Capacity %1\n").arg(p0.PVT.Rso.size()); out << QString("PVT.Bo: Capacity %1\n").arg(p0.PVT.Bo.size()); out << QString("PVT.Co: Capacity %1\n").arg(p0.PVT.Co.size()); out << QString("PVT.miuo: Capacity %1\n").arg(p0.PVT.miuo.size()); out << QString("PVT.rouo: Capacity %1\n").arg(p0.PVT.rouo.size()); out << QString("PVT.Rv: Capacity %1\n").arg(p0.PVT.Rv.size()); out << QString("PVT.Bg: Capacity %1\n").arg(p0.PVT.Bg.size()); out << QString("PVT.Cg: Capacity %1\n").arg(p0.PVT.Cg.size()); out << QString("PVT.miug: Capacity %1\n").arg(p0.PVT.miug.size()); out << QString("PVT.roug: Capacity %1\n").arg(p0.PVT.roug.size()); out << QString("PVT.Z: Capacity %1\n").arg(p0.PVT.Z.size()); out << QString("PVT.Rsw: Capacity %1\n").arg(p0.PVT.Rsw.size()); out << QString("PVT.Bw: Capacity %1\n").arg(p0.PVT.Bw.size()); out << QString("PVT.Cw: Capacity %1\n").arg(p0.PVT.Cw.size()); out << QString("PVT.miuw: Capacity %1\n").arg(p0.PVT.miuw.size()); out << QString("PVT.rouw: Capacity %1\n").arg(p0.PVT.rouw.size()); out << QString("PVT.V: Capacity %1\n").arg(p0.PVT.V.size()); out << QString("PVT.k_kinitial: Capacity %1\n").arg(p0.PVT.k_kinitial.size()); out << QString("PVT.Cf_Cfinitial: Capacity %1\n").arg(p0.PVT.Cf_Cfinitial.size()); out << QString("PVT.So: Capacity %1\n").arg(p0.PVT.So.size()); out << QString("PVT.Kro: Capacity %1\n").arg(p0.PVT.Kro.size()); out << QString("PVT.Sg: Capacity %1\n").arg(p0.PVT.Sg.size()); out << QString("PVT.Krg: Capacity %1\n").arg(p0.PVT.Krg.size()); out << QString("PVT.Sw: Capacity %1\n").arg(p0.PVT.Sw.size()); out << QString("PVT.Krw: Capacity %1\n").arg(p0.PVT.Krw.size()); // 7. Base Data out << "\n--- Base Data ---\n"; out << QString("Base.Pi: %1\n").arg(p0.Base.Pi); out << QString("Base.Cti: %1\n").arg(p0.Base.Cti); out << QString("Base.Cf: %1\n").arg(p0.Base.Cf); out << QString("Base.Soi: %1\n").arg(p0.Base.Soi); out << QString("Base.Sgi: %1\n").arg(p0.Base.Sgi); out << QString("Base.Swi: %1\n").arg(p0.Base.Swi); out << QString("Base.k: Capacity %1\n").arg(p0.Base.k.size()); out << QString("Base.phi: Capacity %1\n").arg(p0.Base.phi.size()); out << QString("Base.h: Capacity %1\n").arg(p0.Base.h.size()); out << QString("Base.d: Capacity %1\n").arg(p0.Base.d); out << QString("Base.dt_Max: Capacity %1\n").arg(p0.Base.dt_Max); out << QString("Base.dt_Min: Capacity %1\n").arg(p0.Base.dt_Min); out << "\n--- End of HX_NWTM_MODEL_INPUT Content ---\n"; file.close(); qDebug() << QString("Successfully saved HX_NWTM_MODEL_INPUT content to: %1").arg(filePath); return true; }