diff --git a/CMakeLists.txt b/CMakeLists.txt index 07377d6..8570f36 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -10,6 +10,7 @@ OPTION(BUILD_SV_INSTALLER "Build SimVascular Installer" OFF) OPTION(buildPy "Build Python Interface" OFF) OPTION(buildDocs "Build Documentation" OFF) OPTION(ENABLE_UNIT_TEST "Enable unit tests" ON) +OPTION(ENABLE_INTERFACE_LIBRARY "Build 1D interface library for 3D-1D coupling" ON) SET(sparseSolverType "skyline" CACHE STRING "Use Sparse Solver") SET_PROPERTY(CACHE sparseSolverType PROPERTY STRINGS skyline superlu csparse) @@ -76,7 +77,7 @@ include(FetchContent) FetchContent_Declare( nlohmann_json GIT_REPOSITORY https://github.com/nlohmann/json.git - GIT_TAG v3.11.2 # You can specify a version or use the latest tag + GIT_TAG v3.12.0 # You can specify a version or use the latest tag ) FetchContent_MakeAvailable(nlohmann_json) # Global include for nlohmann headers @@ -190,3 +191,9 @@ if(ENABLE_UNIT_TEST) add_test(NAME UnitTests COMMAND UnitTestsExecutable WORKING_DIRECTORY ${TESTBIN_DIRECTORY}) endif() + +# BUILD INTERFACE LIBRARY FOR 3D-1D COUPLING +IF(ENABLE_INTERFACE_LIBRARY) + ADD_SUBDIRECTORY("${CMAKE_SOURCE_DIR}/Code/Source/interface") + message(STATUS "Building interface library for 3D-1D coupling") +ENDIF() diff --git a/Code/Source/cvOneDBFSolver.cxx b/Code/Source/cvOneDBFSolver.cxx index 9fbc2c6..8e30b70 100644 --- a/Code/Source/cvOneDBFSolver.cxx +++ b/Code/Source/cvOneDBFSolver.cxx @@ -997,6 +997,376 @@ void cvOneDBFSolver::postprocess_VTK_XML3D_MULTIPLEFILES(){ printf("Results Exported to VTK.\n"); } +// =============================================== +// WRITE 3D XML VTK RESULTS FOR SINGLE TIME STEP +// =============================================== +void cvOneDBFSolver::postprocess_VTK_XML3D_SingleTimeStep( + int timeStep, + cvOneDFEAVector* solution_ptr){ + + // Set Constant Number of Subdivisions on the vessel circumference + int circSubdiv = 20; + + cvOneDSegment* currSeg = NULL; + cvOneDMaterial* curMat = NULL; + cvOneDNode* currNode = NULL; + + // DEFINE INCIDENCE + std::vector segInlets(model->getNumberOfSegments()); + std::vector segOutlets(model->getNumberOfSegments()); + for(int loopSegment=0;loopSegmentgetNumberOfSegments();loopSegment++){ + segInlets[loopSegment] = -1; + segOutlets[loopSegment] = -1; + } + + // FORM INCIDENCE AND STORE COORDS + cvDoubleMat nodeList; + cvDoubleVec temp; + long* segNodes; + int inletNodeID = 0; + int outletNodeID = 0; + if(model->getNumberOfNodes() > 0){ + for(int loopNode = 0; loopNode < model->getNumberOfNodes(); loopNode++){ + temp.clear(); + currNode = model->getNode(loopNode); + temp.push_back(currNode->x); + temp.push_back(currNode->y); + temp.push_back(currNode->z); + nodeList.push_back(temp); + } + for(int loopSegment = 0; loopSegment < model->getNumberOfSegments(); loopSegment++){ + currSeg = model->getSegment(loopSegment); + segNodes = currSeg->getInOutJoints(); + inletNodeID = segNodes[0]; + outletNodeID = segNodes[1]; + segInlets[loopSegment] = inletNodeID; + segOutlets[loopSegment] = outletNodeID; + } + }else{ + currSeg = model->getSegment(0); + temp.clear(); + temp.push_back(0.0); + temp.push_back(0.0); + temp.push_back(0.0); + nodeList.push_back(temp); + temp.clear(); + temp.push_back(currSeg->getSegmentLength()); + temp.push_back(0.0); + temp.push_back(0.0); + nodeList.push_back(temp); + segInlets[0] = 0; + segOutlets[0] = 1; + } + + // Check inlets/outlets + for(int loopSegment=0;loopSegmentgetNumberOfSegments();loopSegment++){ + if(segInlets[loopSegment] == -1){ + printf("ERROR: INLET FOR SEGMENT %d\n",loopSegment); + } + if(segOutlets[loopSegment] == -1){ + printf("ERROR: OUTLET FOR SEGMENT %d\n",loopSegment); + } + } + + // Get solution data from solution_ptr + double* solution_data = solution_ptr->GetEntries(); + + // Set and open VTK file for current time step + string fileName = model->getModelName(); + char timeString[512]; + sprintf(timeString, "_%05d", timeStep); + fileName = fileName + string(timeString) + ".vtp"; + FILE* vtkFile; + vtkFile = fopen(fileName.c_str(),"w"); + + if(vtkFile == NULL){ + cout << "ERROR: Could not open file " << fileName << endl; + return; + } + + // Write VTK XML Header + fprintf(vtkFile,"\n"); + fprintf(vtkFile,"\n"); + fprintf(vtkFile,"\n"); + + // Init Segment Offset + long segOffset = 0; + + // LOOP OVER THE SEGMENTS + for(int loopSegment=0;loopSegmentgetNumberOfSegments();loopSegment++){ + + cvDoubleMat segNodeList; + currSeg = model->getSegment(loopSegment); + + // Compute the total number of points for this segment + int totSegmentSolutions = (currSeg->getNumElements()+1); + int totSegmentPoints = totSegmentSolutions * circSubdiv; + + // Set the range for the solution of this segment + int startOut = segOffset; + int finishOut = segOffset + 2*(totSegmentSolutions); + + // Get Material + curMat = subdomainList[loopSegment]->GetMaterial(); + + // Write Piece Header + fprintf(vtkFile,"\n", + totSegmentPoints, currSeg->getNumElements()); + + // Get inlet and outlet joints + int inletSegJoint = segInlets[loopSegment]; + int outletSegJoint = segOutlets[loopSegment]; + + // Compute Segment Versor + double segVers[3][3]; + segVers[0][0] = nodeList[outletSegJoint][0] - nodeList[inletSegJoint][0]; + segVers[1][0] = nodeList[outletSegJoint][1] - nodeList[inletSegJoint][1]; + segVers[2][0] = nodeList[outletSegJoint][2] - nodeList[inletSegJoint][2]; + double mod = sqrt(segVers[0][0]*segVers[0][0] + segVers[1][0]*segVers[1][0] + segVers[2][0]*segVers[2][0]); + segVers[0][0] /= mod; + segVers[1][0] /= mod; + segVers[2][0] /= mod; + + double lengthByNodes = mod; + double lengthBySegment = currSeg->getSegmentLength(); + + // Compute Segment Local axis system + evalSegmentLocalAxis(segVers); + + // Loop on the number of elements + double currCentre[3] = {0.0}; + double currIniArea = 0.0; + double currIniRad = 0.0; + cvDoubleVec tmp; + + for(int loopEl=0;loopElgetNumElements() + 1;loopEl++){ + currCentre[0] = nodeList[inletSegJoint][0] + loopEl*lengthByNodes/double(currSeg->getNumElements())*segVers[0][0]; + currCentre[1] = nodeList[inletSegJoint][1] + loopEl*lengthByNodes/double(currSeg->getNumElements())*segVers[1][0]; + currCentre[2] = nodeList[inletSegJoint][2] + loopEl*lengthByNodes/double(currSeg->getNumElements())*segVers[2][0]; + + currIniArea = currSeg->getInitInletS() + (loopEl/double(currSeg->getNumElements()))*(currSeg->getInitOutletS() - currSeg->getInitInletS()); + currIniRad = sqrt(currIniArea/M_PI); + + // Loop on the subdivisions + for(int loopSubdiv=0;loopSubdiv\n"); + fprintf(vtkFile,"\n"); + for(int loopA=0;loopA\n"); + fprintf(vtkFile,"\n"); + + // Write Strip Incidence and offset + fprintf(vtkFile,"\n"); + fprintf(vtkFile,"\n"); + for(int loopA=0;loopAgetNumElements();loopA++){ + for(int loopB=0;loopB\n"); + fprintf(vtkFile,"\n"); + for(int loopA=0;loopAgetNumElements();loopA++){ + fprintf(vtkFile,"%d ",(loopA+1)*(circSubdiv*2+2)); + } + fprintf(vtkFile,"\n"); + fprintf(vtkFile,"\n"); + fprintf(vtkFile,"\n"); + + // PRINT OUTPUTS + fprintf(vtkFile,"\n"); + + // PRINT FLOW RATES + fprintf(vtkFile,"\n"); + for(int j=startOut+1;j\n"); + + // PRINT AREA + fprintf(vtkFile,"\n"); + for(int j=startOut;j\n"); + + // PRINT RADIAL DISPLACEMENTS AS VECTORS + fprintf(vtkFile,"\n"); + for(int j=startOut;jgetInitInletS() + (((j-startOut)/2)/double(currSeg->getNumElements()))*(currSeg->getInitOutletS() - currSeg->getInitInletS()); + double newArea = solution_data[j]; + double radDisp = sqrt(newArea/M_PI) - sqrt(iniArea/M_PI); + for(int k=0;k\n"); + + // PRINT PRESSURE IN MMHG + fprintf(vtkFile,"\n"); + double segLength = currSeg->getSegmentLength(); + int section = 0; + for(int j=startOut;jgetNumElements())*segLength; + for(int k=0;kGetPressure(solution_data[j],z)*baryeTommHg); + } + fprintf(vtkFile,"\n"); + section++; + } + fprintf(vtkFile,"\n"); + + // PRINT REYNOLDS NUMBER + fprintf(vtkFile,"\n"); + for(int j=startOut;jGetDensity()/curMat->GetDynamicViscosity()*flo/sqrt(area)*sqrt(4.0/M_PI); + for(int k=0;k\n"); + + // PRINT WSS + fprintf(vtkFile,"\n"); + for(int j=startOut;jGetDynamicViscosity()*flo/(M_PI*radius*radius*radius); + for(int k=0;k\n"); + + // Close Pointdata and Piece + fprintf(vtkFile,"\n"); + fprintf(vtkFile,"\n"); + + // Increment Segment Offset + segOffset += 2*(totSegmentSolutions); + + } // End Segment Loop + + // Close VTK file + fprintf(vtkFile,"\n"); + fprintf(vtkFile,"\n"); + fclose(vtkFile); + +} + +// ===================================================== +// WRITE TEXT RESULTS FOR A SINGLE TIME STEP (COUPLING) +// ===================================================== +void cvOneDBFSolver::postprocess_Text_SingleTimeStep(int timeStep, + cvOneDFEAVector* solution_ptr){ + double* solution_data = solution_ptr->GetEntries(); + int elCount = 0; + + for(int fileIter = 0; fileIter < model->getNumberOfSegments(); fileIter++){ + cvOneDSegment* curSeg = model->getSegment(fileIter); + cvOneDMaterial* curMat = subdomainList[fileIter]->GetMaterial(); + long numEls = curSeg->getNumElements(); + double segLength = curSeg->getSegmentLength(); + long startOut = elCount; + long finishOut = elCount + ((numEls + 1) * 2); + + string base = string(model->getModelName()) + string(curSeg->getSegmentName()); + + string flowFile = base + "_flow.dat"; + string areaFile = base + "_area.dat"; + string pressureFile = base + "_pressure.dat"; + string reFile = base + "_Re.dat"; + string wssFile = base + "_wss.dat"; + + ofstream flow, area, pressure, reynolds, wss; + flow.open(flowFile.c_str(), ios::app); + area.open(areaFile.c_str(), ios::app); + pressure.open(pressureFile.c_str(), ios::app); + reynolds.open(reFile.c_str(), ios::app); + wss.open(wssFile.c_str(), ios::app); + + if(!flow.is_open() || !area.is_open() || !pressure.is_open() || + !reynolds.is_open() || !wss.is_open()){ + printf("ERROR: Could not open txt output files for segment %s at time step %d\n", + curSeg->getSegmentName(), timeStep); + elCount += 2 * (numEls + 1); + continue; + } + + flow.precision(OUTPUT_PRECISION); + area.precision(OUTPUT_PRECISION); + pressure.precision(OUTPUT_PRECISION); + reynolds.precision(OUTPUT_PRECISION); + wss.precision(OUTPUT_PRECISION); + // Output flow for each node (one value per node, one row) + for(int j = startOut + 1; j < finishOut; j += 2){ + flow << solution_data[j] << " "; + } + flow << endl; + // Output area, pressure, Reynolds, WSS for each node (one value per node, one row) + int section = 0; + for(int j = startOut; j < finishOut; j += 2){ + double z = (section / double(numEls)) * segLength; + double val_area = solution_data[j]; + double val_flow = solution_data[j + 1]; + double r = sqrt(val_area / M_PI); + double Re = curMat->GetDensity() / curMat->GetDynamicViscosity() + * val_flow / sqrt(val_area) * sqrt(4.0 / M_PI); + double val_pressure = curMat->GetPressure(val_area, z); + double wssVal = (4.0 * curMat->GetDynamicViscosity() * val_flow) + / (M_PI * r * r * r); + area << val_area << " "; + pressure << val_pressure << " "; + reynolds << Re << " "; + wss << wssVal << " "; + section++; + } + area << endl; + pressure << endl; + reynolds << endl; + wss << endl; + elCount += 2 * (numEls + 1); + flow.close(); + area.close(); + pressure.close(); + reynolds.close(); + wss.close(); + } +} + + + + + // ==================== // MAIN SOLUTION DRIVER // ==================== @@ -1054,6 +1424,45 @@ void cvOneDBFSolver::Solve(void){ } } +// Sub function for the 3D-1D coupling +// first part of the Solve function +void cvOneDBFSolver::Solve_initi(int& systemsize_, char* coupling_types_){ + char errStr[256]; + // First check to make sure we've set a model pointer + // Prior to solution attempt. + if(model == NULL){ + cvOneDError::setErrorNumber(ErrorTypeScope::BAD_VALUE); + strcpy(errStr,"In BFSolver::Solve(...), No model pointer was set prior to solution attempt. Bailing out to avoid a coredump!"); + cvOneDError::setErrorString(errStr); + cvOneDError::CallErrorHandler(); + exit(0); + } + + // Query the model for information, this is where + // the subdomain and material information get passed. + QuerryModelInformation(); + + if(std::string(coupling_types_) == "DIR"){ + DefineMthModels(); + }else if(std::string(coupling_types_) == "NEU"){ + DefineMthModels_cpl_neu(); + } + //DefineMthModels(); + + CreateGlobalArrays(); + + long id; + id = model -> getNumberOfSegments(); + + int i; + for (i=0; i(previousSolution->GetDimension()); + //std::cout << "systemsize_: " << (int)systemsize_ << std::endl; + +} + void cvOneDBFSolver::DefineInletFlow(double* time, double* flrt, int num){ flowTime = new double[num]; flowRate = new double[num]; @@ -1082,6 +1491,23 @@ void cvOneDBFSolver::DefineMthModels(){ AddOneModel(branchM); } +void cvOneDBFSolver::DefineMthModels_cpl_neu(){ + mathModels.clear(); + + cout << "Subdomain No. "<SetInflowRate(flowTime, flowRate, numFlowPts, flowTime[numFlowPts-1]); + //Period = flowTime[numFlowPts-1]; + + cvOneDMthBranchModel* branchM = new cvOneDMthBranchModel(subdomainList, jointList, outletList); + AddOneModel(segM); + AddOneModel(branchM); +} + void cvOneDBFSolver::AddOneModel(cvOneDMthModelBase* model){ mathModels.push_back(model); } @@ -1210,7 +1636,7 @@ void cvOneDBFSolver::QuerryModelInformation(void) subdomain->SetBoundCoronaryValues(time, p_lv,num); cout<<"CORONARY boundary condition"< SetBoundValue(boundV); } @@ -1478,7 +1904,7 @@ void cvOneDBFSolver::GenerateSolution(void){ // Add increment increment->Clear(); - cvOneDGlobal::solver->Solve(*increment); + cvOneDGlobal::solver->Solve(*increment); // solve for the increment: LHS*increment = RHS currentSolution->Add(*increment); @@ -1555,30 +1981,326 @@ void cvOneDBFSolver::GenerateSolution(void){ // Increment Iteration Number iter++; - }// End while + }// End while - checkMass += mathModels[0]->CheckMassBalance() * deltaTime; - cout << " Time = " << currentTime << ", "; - cout << "Mass = " << checkMass << ", "; - cout << "Tot iters = " << std::to_string(iter) << endl; + checkMass += mathModels[0]->CheckMassBalance() * deltaTime; + cout << " Time = " << currentTime << ", "; + cout << "Mass = " << checkMass << ", "; + cout << "Tot iters = " << std::to_string(iter) << endl; - // Save solution if needed - if(step % stepSize == 0){ - sprintf( String2, "%ld", (unsigned long)step); - title = String1 + String2; - currentSolution->Rename(title.data()); + // Save solution if needed + if(step % stepSize == 0){ + sprintf( String2, "%ld", (unsigned long)step); + title = String1 + String2; + currentSolution->Rename(title.data()); - double * tmp = currentSolution -> GetEntries(); - int j; + double * tmp = currentSolution -> GetEntries(); + int j; - for(j=0;j GetDimension(); j++){ - TotalSolution[q][j] = tmp[j]; + for(j=0;j GetDimension(); j++){ + TotalSolution[q][j] = tmp[j]; + } + q++; } - q++; - } - *previousSolution = *currentSolution; - iter_total += iter; - } // End global loop + *previousSolution = *currentSolution; + iter_total += iter; + } // End global time loop cout << "\nAvgerage number of Newton-Raphson iterations per time step = "<<(double)iter_total / (double)maxStep<<"\n"<< endl; } + +void cvOneDBFSolver::InitializeAllEquations() { + cout << "[InitializeAllEquations] Starting equation initialization..." << endl; + + if(mathModels.empty()) { + throw cvException("ERROR: mathModels is empty. DefineMthModels() must be called first."); + } + + if(!previousSolution || !currentSolution) { + throw cvException("ERROR: previousSolution or currentSolution is not initialized."); + } + + *currentSolution = *previousSolution; + int numMath = mathModels.size(); + cout << "[InitializeAllEquations] math model size: " << static_cast(numMath) << endl; + for(int i = 0; i < numMath; i++) { + mathModels[i]->EquationInitialize(previousSolution, currentSolution); + } + + cout << "[InitializeAllEquations] Equation initialization completed successfully" << endl; +} + +// =============================================== +// CONVERT SOLUTION FROM [AREA, FLOW] TO [FLOW, PRESSURE] +// =============================================== +void cvOneDBFSolver::ConvertSolutionToFlowPressure( + cvOneDFEAVector* solution_ptr, + double* solution_vector){ + + int num_sol = solution_ptr->GetDimension(); + int num_nodes = num_sol / 2; + double* tmp = solution_ptr->GetEntries(); + + cout << "Converting solution from [Area, Flow] to [Flow, Pressure]..." << endl; + + // Loop over segments to handle each material properly + long nodeOffset = 0; + + for(int loopSegment = 0; loopSegment < model->getNumberOfSegments(); loopSegment++){ + + // Get Current Segment + cvOneDSegment* currSeg = model->getSegment(loopSegment); + + // Get Material for this segment + cvOneDMaterial* curMat = subdomainList[loopSegment]->GetMaterial(); + + // Get segment properties + double segLength = currSeg->getSegmentLength(); + int numElements = currSeg->getNumElements(); + int numNodesInSegment = numElements + 1; + + // Convert nodes in this segment + for(int i = 0; i < numNodesInSegment; i++){ + // Global node index in solution arrays + int globalNodeIdx = nodeOffset + i; + + // Get area and flow from solution_ptr + // Format: [area1][flow1][area2][flow2]... + double area = tmp[2 * globalNodeIdx]; + double flow = tmp[2 * globalNodeIdx + 1]; + + // Calculate z coordinate for this node + double z = (i / (double)numElements) * segLength; + + // Get pressure from material's constitutive relation + double pressure = curMat->GetPressure(area, z); + + // Store in solution_vector + // Format: [flow1][pressure1][flow2][pressure2]... + solution_vector[2 * globalNodeIdx] = flow; + solution_vector[2 * globalNodeIdx + 1] = pressure; + } + + // Update offset for next segment + nodeOffset += numNodesInSegment; + } + + cout << "Solution conversion completed." << endl; +} + +// =================================================================== +// GET CURRENT SOLUTION FROM CVONEDFEAVE CTOR +// =================================================================== +void cvOneDBFSolver::GetCurrentSolution(double* solution_data, int size) { + if(solution_data == NULL) { + throw cvException("ERROR: solution_data is NULL in GetCurrentSolution"); + } + + if(currentSolution == NULL) { + throw cvException("ERROR: currentSolution is not initialized in GetCurrentSolution"); + } + + // check size + if(size != currentSolution->GetDimension()) { + throw cvException("ERROR: Solution size mismatch in GetCurrentSolution"); + } + + // copy currentSolution to solution_data + double* current_data = currentSolution->GetEntries(); + for(int i = 0; i < size; i++) { + solution_data[i] = current_data[i]; + } + + cout << "[GetCurrentSolution] Current solution extracted successfully" << endl; +} + + +void cvOneDBFSolver::InitializeSolutionFromVector(const double* solution_data, int size) { + if(solution_data == NULL) { + throw cvException("ERROR: solution_data is NULL in InitializeSolutionFromVector"); + } + + if(previousSolution == NULL || currentSolution == NULL) { + throw cvException("ERROR: previousSolution or currentSolution is not initialized. Call CreateGlobalArrays first."); + } + + if(size != previousSolution->GetDimension()) { + throw cvException("ERROR: Solution size mismatch in InitializeSolutionFromVector"); + } + + for(int i = 0; i < size; i++) { + (*previousSolution)[i] = solution_data[i]; + (*currentSolution)[i] = solution_data[i]; + } + + cout << "[InitializeSolutionFromVector] Solution vectors reset successfully" << endl; +} + + +cvOneDFEAVector* cvOneDBFSolver::SolveSingleTimeStep(double currentTimeInput, double interpolated_bc_val) { + currentTime = currentTimeInput; + + clock_t tstart_iter; + clock_t tend_iter; + + // sanity check + // cout << "[SolveSingleTimeStep] current solution: " << endl; + // for (int i = 0; i < currentSolution->GetDimension(); i++) { + // cout << " [" << static_cast(i) << "] = " << (*currentSolution)[i] << endl; + // } + // cout << "[SolveSingleTimeStep] previous solution: " << endl; + // for (int i = 0; i < previousSolution->GetDimension(); i++) { + // cout << " [" << static_cast(i) << "] = " << (*previousSolution)[i] << endl; + // } + + + // Update time-dependent models + int numMath = mathModels.size(); + for(int i = 0; i < numMath; i++) { + mathModels[i]->TimeUpdate(currentTime, deltaTime); + } + + // Initialize Newton-Raphson iteration + int iter = 0; + double normf = 1.0; + double norms = 1.0; + + // Advance current time + currentTime += deltaTime; + + // Newton-Raphson Iterations + while(true) { + tstart_iter = clock(); + + increment->Clear(); + + // Form Newton system + for(int i = 0; i < numMath; i++) { + mathModels[i]->FormNewton(lhs, rhs); + } + + // sanity check + // cout << "[SolveSingleTimeStep] current solution before applying boundary conditions: " << endl; + // for (int i = 0; i < currentSolution->GetDimension(); i++) { + // cout << " [" << static_cast(i) << "] = " << (*currentSolution)[i] << endl; + // } + + // Apply boundary conditions + mathModels[0]->ApplyBoundaryConditions(); + + // Calculate residual norms + if(jointList.size() != 0) { + normf = rhs->Norm(L2_norm, 1, 2, jointList[0]->GetGlobal1stLagNodeID()); + norms = rhs->Norm(L2_norm, 0, 2, jointList[0]->GetGlobal1stLagNodeID()); + } else { + normf = rhs->Norm(L2_norm, 1, 2); + norms = rhs->Norm(L2_norm, 0, 2); + } + + if (std::isnan(norms) || std::isnan(normf)) { + throw cvException("Calculated a NaN for the residual."); + } + + // Check convergence (Newton iteration) + if(iter > 0 && normf < convCriteria && norms < convCriteria) { + cout << " iter: " << std::to_string(iter) << " normf: " << std::to_string(normf) << " norms: " << std::to_string(norms) << endl; + break; + } + + // sanity check + // cout << "[SolveSingleTimeStep] current solution before add increment: " << endl; + // for (int i = 0; i < currentSolution->GetDimension(); i++) { + // cout << " [" << static_cast(i) << "] = " << (*currentSolution)[i] << endl; + // } + + // clear before solve increment + increment->Clear(); + + // Solve linear system: get increment = LHS^-1 * RHS + cvOneDGlobal::solver->Solve(*increment); + currentSolution->Add(*increment); + + // sanity check + // cout << "[SolveSingleTimeStep] current solution after add increment: " << endl; + // for (int i = 0; i < currentSolution->GetDimension(); i++) { + // cout << " [" << static_cast(i) << "] = " << (*currentSolution)[i] << endl; + // } + + // If the area goes less than zero, it tells in which segment the error occurs. + int negArea=0; + if(jointList.size() != 0){ + for (long i= 0; i< jointList[0]->GetGlobal1stLagNodeID();i+=2){ + long elCount = 0; + int fileIter = 0; + //check if area <0 or =nan + if (currentSolution->Get(i) < 0.0 || (currentSolution->Get(i) != currentSolution->Get(i))){ + negArea=1; + while (fileIter < model -> getNumberOfSegments()){ + cvOneDSegment *curSeg = model -> getSegment(fileIter); + long numEls = curSeg -> getNumElements(); + long startOut = elCount; + long finishOut = elCount + ((numEls+1)*2); + char *modelname; + char *segname; + if (startOut <= i && i <= finishOut) { + modelname = model-> getModelName(); + segname = curSeg -> getSegmentName(); + std::string msg = "ERROR: The area of segment '" + std::string(segname) + "' is negative."; + throw cvException(msg.c_str()); + } + elCount += 2*(numEls+1); + fileIter++; + } + } + } + } + + if(negArea==1) { + postprocess_Text(); + assert(0); + } + + // Check for negative areas + if(jointList.size() != 0) { + currentSolution->CheckPositive(0, 2, jointList[0]->GetGlobal1stLagNodeID()); + } else { + currentSolution->CheckPositive(0, 2, currentSolution->GetDimension()); + } + + // Set boundary conditions + mathModels[0]->SetBoundaryConditions_coupled(interpolated_bc_val); + + // sanity check + // cout << "[SolveSingleTimeStep] current solution after apply BC: " << endl; + // for (int i = 0; i < currentSolution->GetDimension(); i++) { + // cout << " [" << static_cast(i) << "] = " << (*currentSolution)[i] << endl; + // } + + + tend_iter = clock(); + + cout << " iter: " << std::to_string(iter) << " normf: " << std::to_string(normf) << " norms: " << std::to_string(norms) << endl; + + if(iter > MAX_NONLINEAR_ITERATIONS) { + cout << "Error: Newton not converged, exceed max iterations" << endl; + break; + } + + iter++; + } + // cout << "Exit Newton loop" << endl; + + // Update previous solution for next time step + *previousSolution = *currentSolution; + + return currentSolution; +} + +void cvOneDBFSolver::extractCplBC(double* solution_data, double& CplValue, char* coupling_types){ + // + mathModels[0]->extractCplBC_model(solution_data, CplValue, coupling_types); +} + +void cvOneDBFSolver::extractCplDOF(int& cpldof, char* coupling_types){ + mathModels[0]->extractCpldof(cpldof, coupling_types); +} \ No newline at end of file diff --git a/Code/Source/cvOneDBFSolver.h b/Code/Source/cvOneDBFSolver.h index 6c30fe9..3beb6ee 100644 --- a/Code/Source/cvOneDBFSolver.h +++ b/Code/Source/cvOneDBFSolver.h @@ -74,6 +74,21 @@ class cvOneDBFSolver{ // Solve the blood flow problem static void Solve(void); + // for 3D-1D coupling + static void Solve_initi(int& systemsize_, char* coupling_types_); + static void InitializeAllEquations(); + static void InitializeSolutionFromVector(const double* solution_data, int size); + static void GetCurrentSolution(double* solution_data, int size); + static cvOneDFEAVector* SolveSingleTimeStep(double currentTimeInput, double interpolated_bc_val); + static void ConvertSolutionToFlowPressure(cvOneDFEAVector* solution_ptr, + double* solution_vector); + static void postprocess_VTK_XML3D_SingleTimeStep(int timeStep, + cvOneDFEAVector* solution_ptr); + static void postprocess_Text_SingleTimeStep(int timeStep, + cvOneDFEAVector* solution_ptr); + static void extractCplBC(double* solution_data, double& CplValue, char* coupling_types); + static void extractCplDOF(int& cpldof, char* coupling_types); + // Get the solution; static double GetSolution(int i, int j){return TotalSolution[i][j];}//IV 082103 @@ -97,6 +112,19 @@ class cvOneDBFSolver{ // Find Segment index given the ID static int getSegmentIndex(int segID); + // for debugging + static double GetPreviousSolution(int index) { + return (previousSolution != nullptr) ? previousSolution->Get(index) : 0.0; + } + static double GetCurrentSolution(int index) { + return (currentSolution != nullptr) ? currentSolution->Get(index) : 0.0; + } + static int GetSubdomainListSize() { return static_cast(subdomainList.size()); } + static int GetMathModelsSize() { return static_cast(mathModels.size()); } + static double* GetFlowTime() { return flowTime; } + static double* GetFlowRate() { return flowRate; } + static long GetNumFlowPts() { return numFlowPts; } + private: // Query Model, Allocate Memory, Set Initial Conditions @@ -108,6 +136,7 @@ class cvOneDBFSolver{ static void GenerateSolution(void); //create MthSegmentModel and MthBranchModel if exists. Also specify inflow profile static void DefineMthModels(void); + static void DefineMthModels_cpl_neu(void); static void AddOneModel(cvOneDMthModelBase* model); static bool wasSet; diff --git a/Code/Source/cvOneDEnums.h b/Code/Source/cvOneDEnums.h index cb112b3..88ab9b2 100644 --- a/Code/Source/cvOneDEnums.h +++ b/Code/Source/cvOneDEnums.h @@ -111,6 +111,7 @@ struct BoundCondTypeScope{ PRESSURE_WAVE = 5, // pressure wave in time RCR = 6, // RCR or Windkessel model, without assuming periodicity CORONARY = 9, // Coronary BC, Jongmin Seo & Hyunjin Kim + COUPLED = 10, // Coupled BC, for 1D-3D couplig, Taeouk Kim NOBOUND }; }; diff --git a/Code/Source/cvOneDModelManager.cxx b/Code/Source/cvOneDModelManager.cxx index e495d69..10d668e 100644 --- a/Code/Source/cvOneDModelManager.cxx +++ b/Code/Source/cvOneDModelManager.cxx @@ -90,6 +90,10 @@ int cvOneDModelManager::CreateSegment(char *segName,long segID, double segLen boundT = BoundCondTypeScope::RCR; }else if(!strcmp(boundType, "CORONARY")){ boundT = BoundCondTypeScope::CORONARY; + }else if(!strcmp(boundType, "CORONARY")){ + boundT = BoundCondTypeScope::CORONARY; + }else if(!strcmp(boundType, "COUPLED")){ + boundT = BoundCondTypeScope::COUPLED; }else{ return CV_ERROR; } @@ -152,16 +156,14 @@ int cvOneDModelManager::CreateSegment(char *segName,long segID, double segLen case BoundCondTypeScope::RCR: seg->setBoundRCRValue(value,num); break; - case BoundCondTypeScope::RESISTANCE: seg->setBoundRCRValue(value,num);//using setBoundRCRValue to set resistance and Pd. break; - case BoundCondTypeScope::CORONARY: seg->SetBoundCoronaryValues(value,time,num); break; - default: + default: // NOBOUND or PRESSURE or FLOW. Now COUPLED is also here seg -> setBoundValue(value[0]); break; } diff --git a/Code/Source/cvOneDMthModelBase.cxx b/Code/Source/cvOneDMthModelBase.cxx index b257cf4..8b7f8e3 100644 --- a/Code/Source/cvOneDMthModelBase.cxx +++ b/Code/Source/cvOneDMthModelBase.cxx @@ -190,6 +190,149 @@ void cvOneDMthModelBase::SetBoundaryConditions(){ }// end for loop } +void cvOneDMthModelBase::SetBoundaryConditions_coupled(double interpolated_value){ + // This is to be called after the solution has been updated in the nonlinear + // loop. It overwrites the components of the current solution so that + // Dirichlet boundary conditions are always satisfied. + + long eqNumbers[2]; // two degress of freedom per node + cvOneDSubdomain* sub; + double InitialPressure; + double currP, currS, resistance; + + // Set up Inlet Dirichlet boundary condition (the default is flow rate) + GetNodalEquationNumbers( 0, eqNumbers, 0); + sub= subdomainList[0]; + switch(cvOneDBFSolver::inletBCtype){ + case BoundCondTypeScope::FLOW: + (*currSolution)[eqNumbers[1]] = GetFlowRate(); + break; + + case BoundCondTypeScope::COUPLED:// Neumann coupling + (*currSolution)[eqNumbers[1]] = interpolated_value; + break; + + case BoundCondTypeScope::PRESSURE_WAVE: + (*currSolution)[eqNumbers[0]] = sub->GetMaterial()->GetArea(GetFlowRate(),0); + InitialPressure = flrt[0]; + break; + default: + + break; + } + + for (vector::iterator it=outletList.begin(); it!=outletList.end(); it++){ + GetNodalEquationNumbers(subdomainList[*it]->GetNumberOfNodes() - 1, eqNumbers, *it); + sub = subdomainList[*it]; + + // Speciafically for RCR BC if treated as essential BC // added IV 051303 + /* double alphaRCR, Rp, Rd, Cap, prevP; + //double InitialQ; + //double MemoI, MemoK, dMemoIdP, dMemoKdP; + //double lhs_QQ, lhs_QS, rhs_Q; + double MemoC; + double z = sub->GetOutletZ(); + */ + switch(sub->GetBoundCondition()){ + case BoundCondTypeScope::PRESSURE: + (*currSolution)[eqNumbers[0]] = sub->GetBoundFlowRate(); + break; + + case BoundCondTypeScope::COUPLED: //Dirichlet coupling + // convert pressure value to area + (*currSolution)[eqNumbers[0]] = sub->GetMaterial()->GetArea(interpolated_value, sub->GetLength()); + break; + + case BoundCondTypeScope::FLOW: + (*currSolution)[eqNumbers[1]] = sub->GetBoundFlowRate(); + break; + case BoundCondTypeScope::RESISTANCE: + if(cvOneDGlobal::CONSERVATION_FORM==0){ + currS = (*currSolution)[eqNumbers[0]]; + currP = sub->GetMaterial()->GetPressure(currS, sub->GetLength()); + resistance = sub->GetBoundResistance(); + (*currSolution)[eqNumbers[1]] = currP/resistance; + } + break; + case BoundCondTypeScope::RESISTANCE_TIME: + if(cvOneDGlobal::CONSERVATION_FORM == 0){ + currS = (*currSolution)[eqNumbers[0]]; + currP = sub->GetMaterial()->GetPressure(currS, sub->GetLength()); + resistance = sub->GetBoundResistance(currentTime); + (*currSolution)[eqNumbers[1]] = currP/resistance; + } + break; + case BoundCondTypeScope::RCR: + break; + default: + break; + }// end switch + + }// end for loop +} + +void cvOneDMthModelBase::extractCplBC_model(double* solution_data, double& CplValue, char* coupling_types){ + // solution_data format: [flow1][pressure1][flow2][pressure2]... for each nodes + long eqNumbers[2]; // two degress of freedom per node + cvOneDSubdomain* sub; + + if(std::string(coupling_types) == "NEU"){ + + if(cvOneDBFSolver::inletBCtype == BoundCondTypeScope::COUPLED) { + GetNodalEquationNumbers( 0, eqNumbers, 0); + CplValue = solution_data[eqNumbers[1]]; // pressure dof for inlet node + + //cout << "CplValue: " << CplValue << endl; + + }else{ + cout << "[ExtractCoupledBCValue] WARNING: NEU coupling but inlet BC is not COUPLED" << endl; + } + + }else if(std::string(coupling_types) == "DIR"){ + for(auto it = outletList.begin(); it != outletList.end(); it++) { + GetNodalEquationNumbers(subdomainList[*it]->GetNumberOfNodes() - 1, eqNumbers, *it); + sub = subdomainList[*it]; + // Check if outlet is coupled + if(sub->GetBoundCondition() == BoundCondTypeScope::COUPLED) { + // Extract flow from solution vector + CplValue = solution_data[eqNumbers[0]]; // flow dof for outlet node + break; // Only extract from first coupled outlet + } + } + } +} + +void cvOneDMthModelBase::extractCpldof(int& Cpldof, char* coupling_types){ + // solution_data format: [flow1][pressure1][flow2][pressure2]... for each nodes + long eqNumbers[2]; // two degress of freedom per node + cvOneDSubdomain* sub; + + if(std::string(coupling_types) == "NEU"){ + + if(cvOneDBFSolver::inletBCtype == BoundCondTypeScope::COUPLED) { + GetNodalEquationNumbers( 0, eqNumbers, 0); + Cpldof = eqNumbers[1]; // pressure dof for inlet node + + //cout << "CplValue: " << CplValue << endl; + + }else{ + cout << "[ExtractCoupledBCValue] WARNING: NEU coupling but inlet BC is not COUPLED" << endl; + } + + }else if(std::string(coupling_types) == "DIR"){ + for(auto it = outletList.begin(); it != outletList.end(); it++) { + GetNodalEquationNumbers(subdomainList[*it]->GetNumberOfNodes() - 1, eqNumbers, *it); + sub = subdomainList[*it]; + // Check if outlet is coupled + if(sub->GetBoundCondition() == BoundCondTypeScope::COUPLED) { + // Extract flow from solution vector + Cpldof = eqNumbers[0]; // flow dof for outlet node + break; // Only extract from first coupled outlet + } + } + } +} + // Eval Mass Balance double cvOneDMthModelBase::CheckMassBalance(){ @@ -236,6 +379,9 @@ void cvOneDMthModelBase::ApplyBoundaryConditions(){ if(cvOneDBFSolver::inletBCtype == BoundCondTypeScope::FLOW){ GetNodalEquationNumbers(0, eqNumbers, 0); cvOneDGlobal::solver->SetSolution(eqNumbers[1], value); + }else if (cvOneDBFSolver::inletBCtype == BoundCondTypeScope::COUPLED){// Neumann coupling for inlet flow BC + GetNodalEquationNumbers(0, eqNumbers, 0); + cvOneDGlobal::solver->SetSolution(eqNumbers[1], value); }else if (cvOneDBFSolver::inletBCtype == BoundCondTypeScope::PRESSURE_WAVE){ GetNodalEquationNumbers(0, eqNumbers, 0); cvOneDGlobal::solver->SetSolution(eqNumbers[0], value); @@ -260,6 +406,7 @@ void cvOneDMthModelBase::ApplyBoundaryConditions(){ switch(sub->GetBoundCondition()){ case BoundCondTypeScope::PRESSURE: case BoundCondTypeScope::PRESSURE_WAVE: + case BoundCondTypeScope::COUPLED: // Dirichlet coupling for outlet pressure BC cvOneDGlobal::solver->SetSolution( eqNumbers[0], value); break; case BoundCondTypeScope::FLOW: @@ -319,6 +466,9 @@ void cvOneDMthModelBase::ApplyBoundaryConditions(){ if(cvOneDBFSolver::inletBCtype == BoundCondTypeScope::FLOW){ GetNodalEquationNumbers( 0, eqNumbers, 0); cvOneDGlobal::solver->SetSolution( eqNumbers[1], value); + }else if (cvOneDBFSolver::inletBCtype == BoundCondTypeScope::COUPLED){ // Neumann coupling for inlet flow BC + GetNodalEquationNumbers(0, eqNumbers, 0); + cvOneDGlobal::solver->SetSolution(eqNumbers[1], value); }else if (cvOneDBFSolver::inletBCtype == BoundCondTypeScope::PRESSURE_WAVE){ GetNodalEquationNumbers( 0, eqNumbers, 0); cvOneDGlobal::solver->SetSolution( eqNumbers[0], value); @@ -388,6 +538,7 @@ void cvOneDMthModelBase::ApplyBoundaryConditions(){ // so same treatment as regular Essential BC like in Brooke's case BoundCondTypeScope::PRESSURE: case BoundCondTypeScope::PRESSURE_WAVE: + case BoundCondTypeScope::COUPLED: // Dirichlet coupling for outlet BC cvOneDGlobal::solver->SetSolution( eqNumbers[0], value); break; case BoundCondTypeScope::FLOW: @@ -628,3 +779,42 @@ double cvOneDMthModelBase::GetFlowRate(){ } +// for Neumann coupled BC, the inlet flow rate is given by the 3D solver as a function of time +double cvOneDMthModelBase::GetCoupledInletFlowRate(const std::vector& params) { + // Expected format: + // params = [2, t1, t2, Q1, Q2] + + if (params.size() != 5) { + throw std::runtime_error( + "GetCoupledInletFlowRate: params size must be 5: [2, t1, t2, Q1, Q2]"); + } + + int num_time_pts = static_cast(params[0]); + if (num_time_pts != 2) { + throw std::runtime_error( + "GetCoupledInletFlowRate: params[0] must be 2"); + } + + double t1 = params[1]; + double t2 = params[2]; + double q1 = params[3]; + double q2 = params[4]; + + if (t2 == t1) { + throw std::runtime_error( + "GetCoupledInletFlowRate: t1 and t2 must be different"); + } + + // Clamp outside the interval + if (currentTime <= t1) { + return q1; + } + if (currentTime >= t2) { + return q2; + } + + // Linear interpolation + double xi = (currentTime - t1) / (t2 - t1); + return q1 + xi * (q2 - q1); +} + diff --git a/Code/Source/cvOneDMthModelBase.h b/Code/Source/cvOneDMthModelBase.h index d996fba..95d8a97 100644 --- a/Code/Source/cvOneDMthModelBase.h +++ b/Code/Source/cvOneDMthModelBase.h @@ -66,6 +66,9 @@ class cvOneDMthModelBase{ // forms minus the global residual vector and an approximation to the global consistent tangent virtual void FormNewton(cvOneDFEAMatrix* lhsMatrix, cvOneDFEAVector* rhsVector) = 0; virtual void SetBoundaryConditions(); + virtual void SetBoundaryConditions_coupled(double interpolated_value); // for 3D-1D coupling + virtual void extractCplBC_model(double* solution_data, double& CplValue, char* coupling_types); // for 3D-1D coupling + virtual void extractCpldof(int& Cpldof, char* coupling_types); // for 3D-1D coupling virtual double CheckMassBalance(); virtual void ApplyBoundaryConditions(); virtual void GetNodalEquationNumbers( long node, long* eqNumbers, long ith); @@ -79,6 +82,7 @@ class cvOneDMthModelBase{ protected: double GetFlowRate(); + double GetCoupledInletFlowRate(const std::vector& params);// for Neumann coupling inlet BC typeOfEquation type; long numberOfEquations; diff --git a/Code/Source/cvOneDMthSegmentModel.cxx b/Code/Source/cvOneDMthSegmentModel.cxx index cf13e87..0f42190 100644 --- a/Code/Source/cvOneDMthSegmentModel.cxx +++ b/Code/Source/cvOneDMthSegmentModel.cxx @@ -80,6 +80,26 @@ void cvOneDMthSegmentModel::FormNewton(cvOneDFEAMatrix* lhsMatrix, cvOneDFEAVect FormElement_FD(element, i, &elementVector, &elementMatrix); rhsVector->Add(elementVector); lhsMatrix->Add(elementMatrix); + + // sanity check + // std::cout << "[FormNewton] Subdomain=" << (int)i + // << " Element=" << (long)element << std::endl; + + // std::cout << " elementVector: "; + // for(int j = 0; j < 4; j++) { + // std::cout << elementVector[j] << " "; + // } + // std::cout << std::endl; + + // double* entries = elementMatrix.GetPointerToEntries(); // ✓ 포인터 얻기 + // std::cout << " elementMatrix (4x4):" << std::endl; + // for(int row = 0; row < 4; row++) { + // std::cout << " "; + // for(int col = 0; col < 4; col++) { + // std::cout << entries[row * 4 + col] << " "; // ✓ 포인터로 접근 + // } + // std::cout << std::endl; + // } } } } @@ -606,7 +626,7 @@ void cvOneDMthSegmentModel::FormElement(long element, } if(bound==BoundCondTypeScope::NOBOUND||bound==BoundCondTypeScope::PRESSURE - ||bound==BoundCondTypeScope::FLOW){ + ||bound==BoundCondTypeScope::FLOW||bound==BoundCondTypeScope::COUPLED){ //Outlet flux term (at z=z_outlet) which is the linearized F-KU if (element == (sub->GetNumberOfElements())-1){ double z = sub->GetOutletZ(); @@ -654,7 +674,7 @@ void cvOneDMthSegmentModel::FormElement(long element, }// end inlet flux if(bound==BoundCondTypeScope::NOBOUND||bound==BoundCondTypeScope::PRESSURE - ||bound==BoundCondTypeScope::FLOW){ + ||bound==BoundCondTypeScope::FLOW||bound==BoundCondTypeScope::COUPLED){ //If no outlet BC or Dirichlet outlet BC, compute the Outlet full flux term (at z=z_outlet) which is the linearized F-KU IV 02-03-03 if (element == (sub->GetNumberOfElements())-1){ double z = sub->GetOutletZ();//checked IV 02-03-03 diff --git a/Code/Source/cvOneDOptions.h b/Code/Source/cvOneDOptions.h index b1a82c9..1f67fd3 100644 --- a/Code/Source/cvOneDOptions.h +++ b/Code/Source/cvOneDOptions.h @@ -116,6 +116,11 @@ struct options{ int useIV; int useStab; + // COUPLING OPTIONS + string couplingStatus; // "ON" or "OFF" + string couplingType; // "DIR" or "NEU" + int couplingSubsteps; // integer value (e.g., 100) + // These are to preserve legacy behavior and are // expected to be eventually migrated into a // post-processing step. diff --git a/Code/Source/cvOneDOptionsLegacySerializer.cxx b/Code/Source/cvOneDOptionsLegacySerializer.cxx index 594c53d..53ac9d8 100644 --- a/Code/Source/cvOneDOptionsLegacySerializer.cxx +++ b/Code/Source/cvOneDOptionsLegacySerializer.cxx @@ -266,6 +266,23 @@ void readOptionsLegacyFormat(string inputFile, options* opts){ throw cvException(string("ERROR: Invalid SOLVEROPTIONS Format. Line " + to_string(lineCount) + "\n").c_str()); } solverOptionsDefined = true; + }else if(upper_string(tokenizedString[0]) == std::string("COUPLING")){ + // printf("Found Coupling Options.\n"); + if(tokenizedString.size() > 4){ + throw cvException(string("ERROR: Too many parameters for COUPLING token. Line " + to_string(lineCount) + "\n").c_str()); + }else if(tokenizedString.size() < 4){ + throw cvException(string("ERROR: Not enough parameters for COUPLING token. Line " + to_string(lineCount) + "\n").c_str()); + } + try{ + // string couplingStatus (ON/OFF) + opts->couplingStatus = tokenizedString[1]; + // string couplingType (DIR/NEU) + opts->couplingType = tokenizedString[2]; + // int couplingSubsteps + opts->couplingSubsteps = atoi(tokenizedString[3].c_str()); + }catch(...){ + throw cvException(string("ERROR: Invalid COUPLING Format. Line " + to_string(lineCount) + "\n").c_str()); + } }else if(upper_string(tokenizedString[0]) == std::string("OUTPUT")){ if(tokenizedString.size() > 3){ throw cvException(string("ERROR: Too many parameters for OUTPUT token. Line " + to_string(lineCount) + "\n").c_str()); diff --git a/Code/Source/interface/CMakeLists.txt b/Code/Source/interface/CMakeLists.txt new file mode 100644 index 0000000..96d4a4f --- /dev/null +++ b/Code/Source/interface/CMakeLists.txt @@ -0,0 +1,81 @@ +# Set the library name +set(lib svoned_interface) + +set(CXXSRCS interface.cpp) +set(HDRS interface.h) + + +if(${sparseSolverType} STREQUAL "skyline") + file(GLOB MAIN_SRC_C "${CMAKE_SOURCE_DIR}/Code/Source/*.c" + "${CMAKE_SOURCE_DIR}/Code/Source/*.cxx") + file(GLOB MAIN_SRC_H "${CMAKE_SOURCE_DIR}/Code/Source/*.h") + +elseif(${sparseSolverType} STREQUAL "superlu") + file(GLOB MAIN_SRC_C "${CMAKE_SOURCE_DIR}/Code/Source/*.c" + "${CMAKE_SOURCE_DIR}/Code/Source/*.cxx" + "${CMAKE_SOURCE_DIR}/Code/Source/sparse/*.cxx" + "${CMAKE_SOURCE_DIR}/Code/Source/sparse/superlu/*.c" + "${CMAKE_SOURCE_DIR}/Code/Source/sparse/superlu/*.cxx") + file(GLOB MAIN_SRC_H "${CMAKE_SOURCE_DIR}/Code/Source/*.h" + "${CMAKE_SOURCE_DIR}/Code/Source/sparse/*.h" + "${CMAKE_SOURCE_DIR}/Code/Source/sparse/superlu/*.h") + +elseif(${sparseSolverType} STREQUAL "csparse") + file(GLOB MAIN_SRC_C "${CMAKE_SOURCE_DIR}/Code/Source/*.c" + "${CMAKE_SOURCE_DIR}/Code/Source/*.cxx" + "${CMAKE_SOURCE_DIR}/Code/Source/sparse/*.cxx" + "${CMAKE_SOURCE_DIR}/Code/Source/sparse/csparse/*.c" + "${CMAKE_SOURCE_DIR}/Code/Source/sparse/csparse/*.cxx") + file(GLOB MAIN_SRC_H "${CMAKE_SOURCE_DIR}/Code/Source/*.h" + "${CMAKE_SOURCE_DIR}/Code/Source/sparse/*.h" + "${CMAKE_SOURCE_DIR}/Code/Source/sparse/csparse/*.h") +endif() + +list(FILTER MAIN_SRC_C EXCLUDE REGEX "main.cxx") +list(FILTER MAIN_SRC_C EXCLUDE REGEX "interface.cpp") + +list(APPEND MAIN_SRC_C "${CMAKE_SOURCE_DIR}/Code/Source/cvOneDSubdomain.cxx") + +# Create shared library +add_library(${lib} SHARED ${CXXSRCS} ${MAIN_SRC_C} ${MAIN_SRC_H}) + +# Include directories +target_include_directories(${lib} PUBLIC + ${CMAKE_SOURCE_DIR}/Code/Source +) + +# Link against required libraries +target_link_libraries(${lib} PUBLIC + nlohmann_json::nlohmann_json +) + +# Link additional libraries for specific sparse solvers +if(${sparseSolverType} STREQUAL "superlu") + find_package(BLAS REQUIRED) + set(SUPERLU_INCLUDE_DIRS "${SUPERLU_DIR}/SRC") + set(SUPERLU_LIBRARY "${SUPERLU_DIR}/lib/libsuperlu_mt_PTHREAD.a") + + target_include_directories(${lib} PRIVATE ${SUPERLU_INCLUDE_DIRS}) + target_link_libraries(${lib} PRIVATE ${BLAS_LIBRARIES} ${SUPERLU_LIBRARY} pthread) +endif() + +# Set output directory +set(_iface_outdir "${CMAKE_BINARY_DIR}/lib") +set_target_properties(${lib} PROPERTIES + RUNTIME_OUTPUT_DIRECTORY "${_iface_outdir}" + LIBRARY_OUTPUT_DIRECTORY "${_iface_outdir}" + ARCHIVE_OUTPUT_DIRECTORY "${_iface_outdir}" + OUTPUT_NAME "svoned_interface" +) + +# Add lib prefix on Windows +if(WIN32) + set_target_properties(${lib} PROPERTIES + PREFIX "lib" + WINDOWS_EXPORT_ALL_SYMBOLS ON + ) +endif() + +message(STATUS "Configuring 1D solver interface library: ${lib}") +message(STATUS "Interface library output directory: ${_iface_outdir}") +message(STATUS "Sparse solver type: ${sparseSolverType}") \ No newline at end of file diff --git a/Code/Source/interface/interface.cpp b/Code/Source/interface/interface.cpp new file mode 100644 index 0000000..e707a66 --- /dev/null +++ b/Code/Source/interface/interface.cpp @@ -0,0 +1,789 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause + +#include "interface.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "cvOneDGlobal.h" +#include "cvOneDModelManager.h" +#include "cvOneDOptions.h" +#include "cvOneDOptionsJsonParser.h" +#include "cvOneDOptionsJsonSerializer.h" +#include "cvOneDOptionsLegacySerializer.h" + +// ===== Resolve operator<< ambiguity ===== +using std::cout; +using std::cerr; +using std::endl; +using std::string; + +////////////////////////////////////////////////////////// +// Helper functions from main.cxx // +////////////////////////////////////////////////////////// + +std::string removeQuotesIfPresent(const std::string& str) { + std::string result = str; + if (!result.empty() && result.front() == '"' && result.back() == '"') { + result = result.substr(1, result.length() - 2); + } + return result; +} + +struct ArgOptions { + std::optional jsonInput = std::nullopt; + std::optional legacyConversionInput = std::nullopt; + std::optional jsonConversionOutput = std::nullopt; +}; + + +cvOneD::options readLegacyOptions(std::string const& inputFile) { + cvOneD::options opts{}; + cvOneD::readOptionsLegacyFormat(inputFile, &opts); + return opts; +} + +std::optional parseArgsAndHandleOptions(int argc, char** argv) { + // Legacy behavior: single input file + if(argc == 2) { + string inputFile{removeQuotesIfPresent(argv[1])}; + return readLegacyOptions(inputFile); + } + + return std::nullopt; +} + +void setOutputGlobals(const cvOneD::options& opts){ + + if(upper_string(opts.outputType) == "TEXT"){ + cvOneDGlobal::outputType = OutputTypeScope::OUTPUT_TEXT; + }else if(upper_string(opts.outputType) == "VTK"){ + cvOneDGlobal::outputType = OutputTypeScope::OUTPUT_VTK; + }else if(upper_string(opts.outputType) == "BOTH"){ + cvOneDGlobal::outputType = OutputTypeScope::OUTPUT_BOTH; + }else{ + throw cvException("ERROR: Invalid OUTPUT Type.\n"); + } + + if(opts.vtkOutputType){ + cvOneDGlobal::vtkOutputType = *opts.vtkOutputType; + if(cvOneDGlobal::vtkOutputType > 1){ + throw cvException("ERROR: Invalid OUTPUT VTK Type.\n"); + } + } + +} + +size_t findJointNodeIndexOrThrow(const auto& jointNodeName, const auto& nodeNames, const auto& jointName){ + // Return the index of "nodeNames" that corresponds to the "jointNodeName" + // or throw a context-specific error. + + auto const iter = std::find(nodeNames.begin(), nodeNames.end(), jointNodeName); + if (iter == nodeNames.end()) { + std::string const errMsg = "ERROR: The node '" + jointNodeName + "' required by joint '" + + jointName + "' was not found in the list of nodes."; + throw cvException(errMsg.c_str()); + } + return std::distance(nodeNames.begin(), iter); +} + +int getDataTableIDFromStringKey(string key){ + bool found = false; + int count = 0; + while((!found)&&(countgetName()); + // Update Counter + if(!found){ + count++; + } + } + if(!found){ + throw cvException(string("ERROR: Cannot find data table entry from string key: " + key + ".\n").c_str()); + return -1; + }else{ + return count; + } +} + +////////////////////////////////////////////////////////// +// Static member data initialization // +////////////////////////////////////////////////////////// + +int OneDSolverInterface::problem_id_count_ = 0; +std::map OneDSolverInterface::interface_list_; + +//----------------------- +// OneDSolverInterface +//----------------------- +/** + * @brief Constructor for OneDSolverInterface. + * + * @param input_file_name The 1D input file name + */ +OneDSolverInterface::OneDSolverInterface(const std::string& input_file_name) + : input_file_name_(input_file_name) { + problem_id_ = problem_id_count_++; + OneDSolverInterface::interface_list_[problem_id_] = this; +} + +/** + * @brief Destructor for OneDSolverInterface. + */ +OneDSolverInterface::~OneDSolverInterface() {} + +////////////////////////////////////////////////////////// +// Callable C interface functions // +////////////////////////////////////////////////////////// + +extern "C" void initialize_1d(const char* input_file, int& problem_id, int& systemSize, + char* coupling_types); + +extern "C" void set_external_step_size_1d(int problem_id, + double external_step_size); + +extern "C" void return_1d_solution(int problem_id, double* solution_1d, int solution_size); + +extern "C" void update_1d_solution(int problem_id, const double* previous_solution_data, int solution_size); + +extern "C" void run_1d_simulation_step_1d(int problem_id, double current_time, int save_time, char* coupling_types, double* params, + double* solution_vector, double& cplBCvalue, char* last_flag, int& error_code); +extern "C" void extract_coupled_dof(int problem_id, int& coupled_dof, char* coupling_types); + +/** + * @brief Initialize the 1D solver interface for 3D-1D coupling. + * + * @param input_file The name of the 1D solver input file (.in format) + * @param problem_id The returned ID used to identify the 1D problem + * @param coupling_types Type of coupling for each surface ("DIR" or "NEU") + */ +void initialize_1d(const char* input_file, int& problem_id, int& systemSize, + char* coupling_types) { + cout << "========== svOneD initialize ==========" << endl; + + try { + std::string input_file_str(input_file); + cout << "[initialize] input_file: " << input_file_str << endl; + + // Create interface object + auto interface = new OneDSolverInterface(input_file_str); + problem_id = interface->problem_id_; + cout << "[initialize] problem_id: " << (int)problem_id << endl; + + // Parse input file and read 1D configuration + char* argv[] = { (char*)"svOneDSolver", const_cast(input_file) }; + auto const simulationOptions = parseArgsAndHandleOptions(2, argv); + + if (simulationOptions) { + cout << "[initialize] Simulation options loaded successfully" << endl; + cout << "[initialize] Model name: " << simulationOptions->modelName << endl; + cout << "[initialize] Time step: " << simulationOptions->timeStep << endl; + cout << "[initialize] Step size: " << simulationOptions->stepSize << endl; + cout << "[initialize] Max steps: " << simulationOptions->maxStep << endl; + + // Store options in interface for later use + interface->model_name_ = simulationOptions->modelName; + interface->time_step_size_ = simulationOptions->stepSize; + interface->max_step_ = simulationOptions->maxStep; + + // Store coupling options in interface + interface->coupling_status_ = simulationOptions->couplingStatus; + interface->coupling_type_ = simulationOptions->couplingType; + interface->coupling_substeps_ = simulationOptions->couplingSubsteps; + + cout << "[initialize] Coupling status: " << interface->coupling_status_ << endl; + cout << "[initialize] Coupling type: " << interface->coupling_type_ << endl; + cout << "[initialize] Coupling substeps: " << static_cast(interface->coupling_substeps_) << endl; + + //// use functions inside runOneDSolver + // Model checking + const cvOneD::options& opts = *simulationOptions; + cvOneD::validateOptions(opts); + cout << "[initialize] validationOptions completed" << endl; + + // Not sure what this function do right now + setOutputGlobals(opts); + cout << "[initialize] setOutputGlobals completed" << endl; + + // Use functions inside createAndRunModel + // only use functions for creating model + cout << "[initialize] creating model: " << opts.modelName << endl; + // Create model manager + cvOneDModelManager* oned = new cvOneDModelManager((char*)opts.modelName.c_str()); + + // Create nodes + printf("Creating Nodes ... \n"); + int totNodes = opts.nodeName.size(); + int nodeError = CV_OK; + for(int loopA = 0; loopA < totNodes; loopA++) { + // Finally Create Joint + nodeError = oned->CreateNode((char*)opts.nodeName[loopA].c_str(), + opts.nodeXcoord[loopA], opts.nodeYcoord[loopA], opts.nodeZcoord[loopA]); + if(nodeError == CV_ERROR) { + throw cvException(string("ERROR: Error Creating NODE " + to_string(loopA) + "\n").c_str()); + } + } + + // Create joints + printf("Creating Joints ... \n"); + int totJoints = opts.jointName.size(); + int jointError = CV_OK; + int* asInlets = nullptr; + int* asOutlets = nullptr; + string currInletName; + string currOutletName; + int jointInletID = 0; + int jointOutletID = 0; + int totJointInlets = 0; + int totJointOutlets = 0; + for(int loopA = 0; loopA < totJoints; loopA++) { + // GET NAMES FOR INLET AND OUTLET + currInletName = opts.jointInletName[loopA]; + currOutletName = opts.jointOutletName[loopA]; + // FIND JOINTINLET INDEX + jointInletID = getListIDWithStringKey(currInletName, opts.jointInletListNames); + if(jointInletID < 0) { + throw cvException(string("ERROR: Cannot Find JOINTINLET for key " + currInletName).c_str()); + } + totJointInlets = opts.jointInletListNumber[jointInletID]; + // FIND JOINTOUTLET INDEX + jointOutletID = getListIDWithStringKey(currOutletName, opts.jointOutletListNames); + if(jointInletID < 0) { + throw cvException(string("ERROR: Cannot Find JOINTOUTLET for key " + currOutletName).c_str()); + } + // GET TOTALS + totJointInlets = opts.jointInletListNumber[jointInletID]; + totJointOutlets = opts.jointOutletListNumber[jointOutletID]; + // ALLOCATE INLETS AND OUTLET LIST + asInlets = nullptr; + asOutlets = nullptr; + if(totJointInlets > 0) { + asInlets = new int[totJointInlets]; + for(int loopB = 0; loopB < totJointInlets; loopB++) { + asInlets[loopB] = opts.jointInletList[jointInletID][loopB]; + } + } + if(totJointOutlets > 0) { + asOutlets = new int[totJointOutlets]; + for(int loopB = 0; loopB < totJointOutlets; loopB++) { + asOutlets[loopB] = opts.jointOutletList[jointOutletID][loopB]; + } + } + + // Find the index of the indicated node. + auto const jointName = opts.jointName.at(loopA); + auto const nodeIndex = findJointNodeIndexOrThrow( + opts.jointNode.at(loopA), opts.nodeName, jointName); + + // Finally Create Joint + jointError = oned->CreateJoint(jointName.c_str(), + opts.nodeXcoord[nodeIndex], opts.nodeYcoord[nodeIndex], opts.nodeZcoord[nodeIndex], + totJointInlets, totJointOutlets, asInlets, asOutlets); + if(jointError == CV_ERROR) { + throw cvException(string("ERROR: Error Creating JOINT " + to_string(loopA) + "\n").c_str()); + } + // Deallocate + delete[] asInlets; + delete[] asOutlets; + asInlets = nullptr; + asOutlets = nullptr; + } + + // Create materials + printf("Creating Materials ... \n"); + int totMaterials = opts.materialName.size(); + int matError = CV_OK; + double doubleParams[3]; + int matID = 0; + string currMatType = "MATERIAL_OLUFSEN"; + int numParams = 0; + for(int loopA = 0; loopA < totMaterials; loopA++) { + if(upper_string(opts.materialType[loopA]) == "OLUFSEN") { + currMatType = "MATERIAL_OLUFSEN"; + numParams = 3; + } else { + currMatType = "MATERIAL_LINEAR"; + numParams = 1; + } + doubleParams[0] = opts.materialParam1[loopA]; + doubleParams[1] = opts.materialParam2[loopA]; + doubleParams[2] = opts.materialParam3[loopA]; + // CREATE MATERIAL + matError = oned->CreateMaterial((char*)opts.materialName[loopA].c_str(), + (char*)currMatType.c_str(), + opts.materialDensity[loopA], + opts.materialViscosity[loopA], + opts.materialExponent[loopA], + opts.materialPRef[loopA], + numParams, doubleParams, + &matID); + if(matError == CV_ERROR) { + throw cvException(string("ERROR: Error Creating MATERIAL " + to_string(loopA) + "\n").c_str()); + } + } + + // Create datatables + printf("Creating Data Tables ... \n"); + int totCurves = opts.dataTableName.size(); + int curveError = CV_OK; + for(int loopA = 0; loopA < totCurves; loopA++) { + curveError = oned->CreateDataTable((char*)opts.dataTableName[loopA].c_str(),(char*)opts.dataTableType[loopA].c_str(), opts.dataTableVals[loopA]); + if(curveError == CV_ERROR) { + throw cvException(string("ERROR: Error Creating DATATABLE " + to_string(loopA) + "\n").c_str()); + } + } + + // Create segments + printf("Creating Segments ... \n"); + int segmentError = CV_OK; + int totalSegments = opts.segmentName.size(); + int curveTotals = 0; + double* curveTime = nullptr; + double* curveValue = nullptr; + string matName; + string curveName; + int currMatID = 0; + int dtIDX = 0; + for(int loopA = 0; loopA < totalSegments; loopA++) { + + // GET MATERIAL + matName = opts.segmentMatName[loopA]; + currMatID = getListIDWithStringKey(matName, opts.materialName); + if(currMatID < 0) { + throw cvException(string("ERROR: Cannot Find Material for key " + matName).c_str()); + } + + // GET CURVE DATA + curveName = opts.segmentDataTableName[loopA]; + + if(upper_string(curveName) != "NONE") { + dtIDX = getDataTableIDFromStringKey(curveName); + curveTotals = cvOneDGlobal::gDataTables[dtIDX]->getSize(); + curveTime = new double[curveTotals]; + curveValue = new double[curveTotals]; + for(int loopA = 0; loopA < curveTotals; loopA++) { + curveTime[loopA] = cvOneDGlobal::gDataTables[dtIDX]->getTime(loopA); + curveValue[loopA] = cvOneDGlobal::gDataTables[dtIDX]->getValues(loopA); + } + } else { + curveTotals = 1; + curveTime = new double[curveTotals]; + curveValue = new double[curveTotals]; + curveTime[0] = 0.0; + curveValue[0] = 0.0; + } + segmentError = oned->CreateSegment((char*)opts.segmentName[loopA].c_str(), + (long)opts.segmentID[loopA], + opts.segmentLength[loopA], + (long)opts.segmentTotEls[loopA], + (long)opts.segmentInNode[loopA], + (long)opts.segmentOutNode[loopA], + opts.segmentInInletArea[loopA], + opts.segmentInOutletArea[loopA], + opts.segmentInFlow[loopA], + currMatID, + (char*)opts.segmentLossType[loopA].c_str(), + opts.segmentBranchAngle[loopA], + opts.segmentUpstreamSegment[loopA], + opts.segmentBranchSegment[loopA], + (char*)opts.segmentBoundType[loopA].c_str(), + curveValue, + curveTime, + curveTotals); + if(segmentError == CV_ERROR) { + throw cvException(string("ERROR: Error Creating SEGMENT " + to_string(loopA) + "\n").c_str()); + } + // Deallocate + delete[] curveTime; + curveTime = nullptr; + delete[] curveValue; + curveValue = nullptr; + }// until this part of the code is the part of the code of createAndUnModel before SOLVE MODEL + + + string inletCurveName = opts.inletDataTableName; + int inletCurveTotals; + double* inletCurveTime; + double* inletCurveValue; + + if (upper_string(inletCurveName) != "NONE") { // Dirichlet coupling + + int inletCurveIDX = getDataTableIDFromStringKey(inletCurveName); + inletCurveTotals = cvOneDGlobal::gDataTables[inletCurveIDX]->getSize(); + + inletCurveTime = new double[inletCurveTotals]; + inletCurveValue = new double[inletCurveTotals]; + + if (std::string(coupling_types) == "DIR") { + printf("Dirichlet coupling. Get inlet flow data\n"); + + for (int loopB = 0; loopB < inletCurveTotals; loopB++) { + inletCurveTime[loopB] = + cvOneDGlobal::gDataTables[inletCurveIDX]->getTime(loopB); + inletCurveValue[loopB] = + cvOneDGlobal::gDataTables[inletCurveIDX]->getValues(loopB); + } + } + } else { // Neumann coupling + inletCurveTotals = 1; + inletCurveTime = new double[1]; + inletCurveValue = new double[1]; + + inletCurveTime[0] = 0.0; + inletCurveValue[0] = 0.0; + } + + + // 여기부터는 SolveModel중 초기화에 해당되는 부분을 추가 + int solveError = CV_OK; + // inlet boundary type이 opts.boundaryType에 저장되어있는데 (.in 파일에서 SOLVEROPTIONS에서 읽은것) + // 아래 코드는 SolveModel에서 그대로 가지고 온건데 왜 여러 다른 boundary type이 있는지 아직 잘 모르겠음. + // TODO: 나중에 Neumann coupling의 경우 이 inlet BC이 coupled 같은게 될것이므로 수정 필요 + // 일딴은 그냥 넣어보자 + BoundCondTypeScope::BoundCondType boundT; + + // set the creation flag to off. + cvOneDGlobal::isCreating = false; + char* boundType_tmp = (char*)opts.boundaryType.c_str(); + // convert char string to boundary condition type + if(!strcmp( boundType_tmp, "NOBOUND")){ + boundT = BoundCondTypeScope::NOBOUND; + printf("Inlet Condition Type: NOBOUND\n"); + }else if(!strcmp( boundType_tmp, "PRESSURE")){ + boundT = BoundCondTypeScope::PRESSURE; + printf("Inlet Condition Type: PRESSURE\n"); + }else if(!strcmp( boundType_tmp, "PRESSURE_WAVE")){ + boundT = BoundCondTypeScope::PRESSURE_WAVE; + printf("Inlet Condition Type: PRESSURE_WAVE\n"); + }else if(!strcmp( boundType_tmp, "FLOW")){ + boundT = BoundCondTypeScope::FLOW; + printf("Inlet Condition Type: FLOW\n"); + }else if(!strcmp( boundType_tmp, "RESISTANCE")){ + boundT = BoundCondTypeScope::RESISTANCE; + printf("Inlet Condition Type: RESISTANCE\n"); + }else if(!strcmp( boundType_tmp, "RESISTANCE_TIME")){ + boundT = BoundCondTypeScope::RESISTANCE_TIME; + printf("Inlet Condition Type: RESISTANCE_TIME\n"); + }else if(!strcmp( boundType_tmp, "RCR")){ + boundT = BoundCondTypeScope::RCR; + printf("Inlet Condition Type: RCR\n"); + }else if(!strcmp( boundType_tmp, "CORONARY")){ + boundT = BoundCondTypeScope::CORONARY; + printf("Inlet Condition Type: CORONARY\n"); + }else if(!strcmp(boundType_tmp, "COUPLED")){ + boundT = BoundCondTypeScope::COUPLED; + printf("Inlet Condition Type: COUPLED\n"); + }else{ + solveError = CV_ERROR; + } + + // Set Solver Options + cvOneDMthSegmentModel::STABILIZATION = opts.useStab; // 1=stabilization, 0=none + cvOneDGlobal::CONSERVATION_FORM = opts.useIV; + cvOneDBFSolver::ASCII = 1; + + // enroll model inside solver of cvOneDBFSolver file + cvOneDBFSolver::SetModelPtr(cvOneDGlobal::gModelList[cvOneDGlobal::currentModel]); + + // We need to get these from the solver + cvOneDBFSolver::SetDeltaTime(opts.timeStep); + cvOneDBFSolver::SetStepSize(opts.stepSize); + cvOneDBFSolver::SetMaxStep(opts.maxStep); + cvOneDBFSolver::SetQuadPoints(opts.quadPoints); + cvOneDBFSolver::SetInletBCType(boundT); + if(std::string(coupling_types) == "DIR"){// Inflow only need for Dirichlet coupling + cvOneDBFSolver::DefineInletFlow(inletCurveTime, inletCurveValue, inletCurveTotals); + } + cvOneDBFSolver::SetConvergenceCriteria(opts.convergenceTolerance); + + cvOneDGlobal::isSolving = true; + // 여기까지가 SolveModel에서 Solve() 함수 이전까지의 부분 + + // 여기부터는 Solve() 함수 부분 중에서도 GenerateSolution 이전부분 + cvOneDBFSolver::Solve_initi(systemSize, coupling_types); // TODO: 지금 이 안에도 커플링위해서 바꿔야하는 함수들 많음 + // output: systemSize, which is the total number of unknowns in the system. #NODE x 2 (flow&area) + // cout << "system size: " << static_cast(systemSize) << endl; // total number of unknows in the system. #NODE x 2 (flow&area) + + // time loop 시작 전에 필요한 초기화 작업들 + // 여기부터는 Solve() 함수 부분 중에서도 GenerateSolution 이전부분 + cout << "[initialize_1d] Model initialized, preparing for time-stepping..." << endl; + + //EquationInitialize() 호출 + try { + cvOneDBFSolver::InitializeAllEquations(); + cout << "[initialize_1d] Equations initialized successfully" << endl; + } catch (const std::exception& e) { + cerr << "[initialize_1d] Error initializing equations: " << e.what() << endl; + throw; + } + + // Time loop initialization + interface->time_step_ = 0; + + cout << "[initialize_1d] 1D model is ready for time-stepping" << endl; + + }else { + cout << "[initialize_1d] WARNING: No simulation options found" << endl; + } + cout << "[initialize_1d] 1D model initialized successfully" << endl; + + } catch (const std::exception& e) { + cerr << "Error in initialize_1d: " << e.what() << endl; + problem_id = -1; + } +} + +/** + * @brief Set the external (3D solver) time step size for the 1D interface. + * + * @param problem_id The ID used to identify the 1D problem. + * @param external_step_size The time step size of the external program. + */ +void set_external_step_size_1d(int problem_id, double external_step_size) { + auto it = OneDSolverInterface::interface_list_.find(problem_id); + if (it == OneDSolverInterface::interface_list_.end()) { + cerr << "Error: problem_id " << (int)problem_id << " not found" << endl; + return; + } + + + auto interface = it->second; + + double oned_step_size = external_step_size/ (double(interface->coupling_substeps_)); + // interface->coupling_substeps_ : number of sub steps that 1D solver will take within one external step. + // For example, if coupling_substeps_ = 2, then 1D solver will take 2 steps within one external step + // so each 1D step size will be external_step_size/2. + + interface->external_step_size_ = oned_step_size;//now 3D and 1D has same step size + + cvOneDBFSolver::SetDeltaTime(oned_step_size);// set deltaTime in solver as oned_step_size + + + cout << "[set_external_step_size_1d] Step size: " << oned_step_size << " s" << endl; +} + +/** + * @brief copy current 1D solution from 1D solver to 3D solver + * + * @param problem_id The ID used to identify the 1D problem. + * @param solution_1d 1D solution of 3D solver + */ +void return_1d_solution(int problem_id, double* solution_1d, int solution_size){ + auto it = OneDSolverInterface::interface_list_.find(problem_id); + if (it == OneDSolverInterface::interface_list_.end()) { + cerr << "[return_1d_solution] Error: problem_id " + << static_cast(problem_id) << " not found" << endl; + return; + } + + auto interface = it->second; + + try { + cout << "[return_1d_solution] ========================================" << endl; + cout << "[return_1d_solution] Extracting solution for problem_id: " + << static_cast(problem_id) << endl; + + // cvOneDBFSolver의 getter 함수 호출 + cvOneDBFSolver::GetCurrentSolution(solution_1d, solution_size); + + cout << "[return_1d_solution] Solution extracted successfully" << endl; + + } catch (const std::exception& e) { + cerr << "[return_1d_solution] Exception caught: " << e.what() << endl; + throw; + } +} + +/** + * @brief Reset the 1D solver solution vectors for coupling with 3D solver. + * + * This function initializes previousSolution and currentSolution with values from + * the previous time step, which is required when the 1D solver is called multiple times + * within a 3D Newton loop. + * + * @param problem_id The ID of the 1D problem + * @param previous_solution_data Array containing the solution from previous time step + * @param solution_size The size of the solution array + */ +void update_1d_solution(int problem_id, const double* previous_solution_data, int solution_size) { + auto it = OneDSolverInterface::interface_list_.find(problem_id); + if (it == OneDSolverInterface::interface_list_.end()) { + cerr << "[reset_1d_solution] Error: problem_id " + << static_cast(problem_id) << " not found" << endl; + return; + } + + auto interface = it->second; + + try { + cout << "[update_1d_solution] ========================================" << endl; + + // call function from cvOneDBFSolver + cvOneDBFSolver::InitializeSolutionFromVector(previous_solution_data, solution_size); + + cout << "[update_1d_solution] Solution vectors reset successfully" << endl; + + } catch (const std::exception& e) { + cerr << "[update_1d_solution] Exception caught: " << e.what() << endl; + throw; + } +} + + +/** + * @brief Run one time step of the 1D simulation. + * + * Performs a single time step of 1D blood flow simulation by calling the + * underlying solver's SolveSingleTimeStep function. Returns the solution + * vector (pressure and flow) after the time step. + * + * @param problem_id The ID used to identify the 1D problem. + * @param current_time Current simulation time (seconds) + * @param save_time Save 1D results every save_time 3D time steps + * @param coupling_types Type of coupling for each surface ("DIR" or "NEU") + * @param params bc data from 3D solver for interpolation [2, t1, t2, value1, value2] + * @param solution_vector Output vector containing pressure and flow values + * Size: 2*system_size (pressure[0..system_size-1], flow[system_size..2*system_size-1]) + * @param error_code Output error code (0 = success, <0 = error) + */ +void run_1d_simulation_step_1d(int problem_id, double current_time, int save_time, char* coupling_types, double* params, + double* solution_vector, double& cplBCvalue, char* last_flag,int& error_code) { + auto it = OneDSolverInterface::interface_list_.find(problem_id); //problem_id에 해당되는 interface 객체를 찾음 + if (it == OneDSolverInterface::interface_list_.end()) { + cerr << "[run_1d_simulation_step_1d] Error: problem_id " + << static_cast(problem_id) << " not found" << endl; + error_code = -1; + return; + } + + auto interface = it->second; // map에서 problem_id에 해당하는 interface 객체를 가져옴(second: value) + error_code = 0; + auto oned_substeps = interface->coupling_substeps_; + + try { + cout << "[run_1d_simulation_step_1d] ========================================" << endl; + cout << "[run_1d_simulation_step_1d] Time step " + << static_cast(interface->time_step_) + << ", current_time = " << current_time << " s" << endl; + + + // Update current time in solver + cvOneDBFSolver::currentTime = current_time; + + // now params[0] is always 2. not used here + cvOneDFEAVector* solution_ptr = nullptr; + double t1 = params[1]; + double t2 = params[2]; + double val1 = params[3]; + double val2 = params[4]; + double alpha = 0.0; + double interpolated_value = 0.0; + + for(int i = 0; i < oned_substeps; i++) { + // get coupled data from 3D solver and do interpolation for current time + // this is same for Dir or Neu coupling + if (current_time <= t1) { + interpolated_value = val1; + } else if (current_time >= t2) { + interpolated_value = val2; + } else { + alpha = (current_time - t1) / (t2 - t1); + interpolated_value = val1 + alpha * (val2 - val1); + } + // cout << "[run_1d_simulation_step_1d] Substep " << (i+1) << "/" << static_cast(oned_substeps) + // << ", interpolated_value = " << interpolated_value << endl; + + // Call the underlying 1D solver to compute one time step + // SolveSingleTimeStep returns a pointer to the solution vector + solution_ptr = cvOneDBFSolver::SolveSingleTimeStep(current_time, interpolated_value); + + current_time += interface->external_step_size_; // update current_time for next 1D substep + // external step size is 1D solver time step size + } + + + // Check if the solver returned a valid solution + if (solution_ptr == nullptr) { + throw std::runtime_error("SolveSingleTimeStep returned null pointer"); + } + + // Step 1: Copy solution_ptr directly to solution_vector [area1][flow1][area2][flow2]... + int solution_size = solution_ptr->GetDimension(); + double* solution_data_tmp = solution_ptr->GetEntries(); + for(int i = 0; i < solution_size; i++) { + solution_vector[i] = solution_data_tmp[i]; + } + + // Step 2: Convert to flow and pressure format separately (for internal use or debugging) + // Create a separate vector for converted solution + // Make solution vector to transfer to 3D solver + // converted_solution format: [flow1][pressure1][flow2][pressure2]... for each nodes + // currnet solution_ptr: [area1][flow1][area2][flow2]... for each nodes + double* converted_solution = new double[solution_size]; + cvOneDBFSolver::ConvertSolutionToFlowPressure(solution_ptr, converted_solution); + // converted_solution now contains [flow1][pressure1][flow2][pressure2]... + // You can use this for other purposes if needed + + // extract coupled BC value for 3D solver from converted_solution + cvOneDBFSolver::extractCplBC(converted_solution, cplBCvalue, coupling_types); + + delete[] converted_solution; + + + // last_flag is "L" when it is at the last iteration + // last_flag is "D: when it is at the middel of 3D newton iteration + if (std::string(last_flag) == "L"){ + // Update interface internal states + interface->time_step_++; // this is for 1D internal use only. substep is not inlcuded. + // this is same time_step with 3D solver + + // print solution as vtk file and txt files + if (interface->time_step_ % save_time == 0) { + cout << "generate vtk file at time step: "<< static_cast(interface->time_step_) << endl; + cvOneDBFSolver::postprocess_VTK_XML3D_SingleTimeStep(interface->time_step_, solution_ptr); + cout << "generate txt files at time step: "<< static_cast(interface->time_step_) << endl; + cvOneDBFSolver::postprocess_Text_SingleTimeStep(interface->time_step_, solution_ptr); + } + + } + + + + cout << "[run_1d_simulation_step_1d] Time step completed" << endl; + + } catch (const std::exception& e) { + cerr << "[run_1d_simulation_step_1d] Error: " << e.what() << endl; + error_code = -1; + } +} + +void extract_coupled_dof(int problem_id, int& coupled_dof, char* coupling_types){ + auto it = OneDSolverInterface::interface_list_.find(problem_id); + if (it == OneDSolverInterface::interface_list_.end()) { + cerr << "[extract_coupled_dof] Error: problem_id " + << static_cast(problem_id) << " not found" << endl; + return; + } + + auto interface = it->second; + + try { + cout << "[extract_coupled_dof] ========================================" << endl; + + // call function from cvOneDBFSolver + cvOneDBFSolver::extractCplDOF(coupled_dof, coupling_types); + + cout << "[extract_coupled_dof] Coupled DOF extracted successfully" << endl; + + } catch (const std::exception& e) { + cerr << "[extract_coupled_dof] Exception caught: " << e.what() << endl; + throw; + } + +} \ No newline at end of file diff --git a/Code/Source/interface/interface.h b/Code/Source/interface/interface.h new file mode 100644 index 0000000..ba64ba8 --- /dev/null +++ b/Code/Source/interface/interface.h @@ -0,0 +1,71 @@ +/** + * @file interface.h + * @brief svOneDSolver callable interface for 3D-1D coupling. + */ + +#include +#include +#include + +/** + * @brief Interface class for calling svOneD from external programs + */ +class OneDSolverInterface { + public: + /** + * @brief Construct a new interface object + * @param input_file_name The 1D JSON file which specifies the model + */ + OneDSolverInterface(const std::string& input_file_name); + + /** + * @brief Destroy the interface object + */ + ~OneDSolverInterface(); + + /** + * @brief Counter for the number of interfaces + */ + static int problem_id_count_; + + /** + * @brief List of interfaces + */ + static std::map interface_list_; + + /** + * @brief ID of current interface + */ + int problem_id_ = 0; + + /** + * @brief 1D input (JSON) file + */ + std::string input_file_name_; + + /** + * @brief Time step size of the external program (3D solver) + * + * This is required for coupling with a 3D solver + */ + double external_step_size_ = 0.001; + + /** + * @brief Current time step + */ + int time_step_ = 0; + + /** + * @brief Current simulation time (seconds) + */ + double current_time_ = 0.0; + + std::string model_name_; + double time_step_size_ = 0.0; + int max_step_ = 0; + + // Coupling options + std::string coupling_status_ = "OFF"; // "ON" or "OFF" + std::string coupling_type_; // "DIR" or "NEU" + int coupling_substeps_ = 10; // number of substeps +}; \ No newline at end of file diff --git a/Code/Source/main.cxx b/Code/Source/main.cxx index 916b636..9b92281 100644 --- a/Code/Source/main.cxx +++ b/Code/Source/main.cxx @@ -292,12 +292,13 @@ void createAndRunModel(const cvOneD::options& opts) { curveValue = nullptr; } - double* vals; - int tot; + double* vals;// not used + int tot;// not used // SOLVE MODEL printf("Solving Model ... \n"); int solveError = CV_OK; + // Get Inlet BC data string inletCurveName = opts.inletDataTableName; int inletCurveIDX = getDataTableIDFromStringKey(inletCurveName); int inletCurveTotals = cvOneDGlobal::gDataTables[inletCurveIDX]->getSize(); @@ -460,10 +461,11 @@ cvOneD::options readLegacyOptions(std::string const& inputFile){ // std::optional parseArgsAndHandleOptions(int argc, char** argv){ - // Legacy behavior + // Legacy behavior. right now, argc==2. + // Legacy input file: "xxx.in" stored in argv[1]. if(argc == 2){ string inputFile{removeQuotesIfPresent(argv[1])}; - return readLegacyOptions(inputFile); + return readLegacyOptions(inputFile); //ends at here } // Default behavior @@ -493,7 +495,7 @@ int main(int argc, char** argv){ try{ - auto const simulationOptions = parseArgsAndHandleOptions(argc,argv); + auto const simulationOptions = parseArgsAndHandleOptions(argc,argv); // read file if(simulationOptions){ // The simulation options were defined so we can run the simulation runOneDSolver(*simulationOptions);