diff --git a/tools/getFromFoam/getFromFoam.py b/tools/getFromFoam/getFromFoam.py deleted file mode 100644 index e69de29..0000000 diff --git a/tools/linAlg4Foam/OpenFoamFunctionObjects/FOL2norm b/tools/linAlg4Foam/OpenFoamFunctionObjects/FOL2norm new file mode 100644 index 0000000..3519c04 --- /dev/null +++ b/tools/linAlg4Foam/OpenFoamFunctionObjects/FOL2norm @@ -0,0 +1,104 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: v2406 | +| \\ / A nd | Website: www.openfoam.com | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ + +// This function object calculates the L2 norm of a field. +// +// Input arguments: +// - path to the snapshot (time folder) +// - field + +L2norm +{ + type coded; + libs (utilityFunctionObjects); + region region0; + + name L2norm; + + codeWrite + #{ + // Read the list of fields in python format (e.g., ['T', 'DT']) + string inputLine; + std::cout << "Enter list of fields: "; + std::getline(std::cin, inputLine); + + // Strip square brackets + if (!inputLine.empty() && inputLine.front() == '[') inputLine.erase(0, 1); + if (!inputLine.empty() && inputLine.back() == ']') inputLine.pop_back(); + + // Split and clean elements + wordList fieldList; + std::stringstream ss(inputLine); + std::string item; + + while (std::getline(ss, item, ',')) { + // Remove single quotes and whitespace + item.erase(std::remove(item.begin(), item.end(), '\''), item.end()); + item.erase(0, item.find_first_not_of(" \t")); + item.erase(item.find_last_not_of(" \t") + 1); + + if (item.find("Region") != std::string::npos) { + item = "../" + item; + } + + fieldList.append(word(item)); + } + + Info << "Parsed fields:" << nl; + forAll(fieldList, i) + { + Info << " - " << fieldList[i] << nl; + } + + // Loop for field + snapshot sampling + while (true) + { + // Define strings + string snapshotPath; + + // Get snapshot path + std::cout << "Enter snapshot path: "; + if (!std::getline(std::cin, snapshotPath) || snapshotPath.empty()) + { + break; + } + + // Construct snapshot path + std::string fieldTime = snapshotPath; + + forAll(fieldList, fieldI) + { + volScalarField field + ( + IOobject + ( + fieldList[fieldI], + fieldTime, + mesh(), + IOobject::MUST_READ, + IOobject::NO_WRITE + ), + mesh() + ); + + volScalarField field2 = magSqr(field); + dimensionedScalar L2sum = Foam::fvc::domainIntegrate(field2); + + scalar L2norm = sqrt(L2sum.value()); + + // Info output + string snapshotID = snapshotPath.substr(snapshotPath.rfind('/', snapshotPath.rfind('/') - 1) + 1); + Info << "L2norm: " << snapshotID << ", " << fieldList[fieldI] << " = " << L2norm << endl; + } + } + + #}; +} + +// ************************************************************************* // + diff --git a/tools/linAlg4Foam/OpenFoamFunctionObjects/FOfieldMinMax b/tools/linAlg4Foam/OpenFoamFunctionObjects/FOfieldMinMax new file mode 100644 index 0000000..0b90ee9 --- /dev/null +++ b/tools/linAlg4Foam/OpenFoamFunctionObjects/FOfieldMinMax @@ -0,0 +1,127 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: v2406 | +| \\ / A nd | Website: www.openfoam.com | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ + +// This function object calculates the min/max values and their positions +// in field. +// +// Input arguments: +// - path to the snapshot (time folder) +// - field +// + +fieldMinMax +{ + type coded; + libs (utilityFunctionObjects); + region region0; + + name fieldMinMax; + + codeWrite + #{ + // Read the list of fields in python format (e.g., ['T', 'DT']) + string inputLine; + std::cout << "Enter list of fields: "; + std::getline(std::cin, inputLine); + + // Strip square brackets + if (!inputLine.empty() && inputLine.front() == '[') inputLine.erase(0, 1); + if (!inputLine.empty() && inputLine.back() == ']') inputLine.pop_back(); + + // Split and clean elements + wordList fieldList; + std::stringstream ss(inputLine); + std::string item; + + while (std::getline(ss, item, ',')) { + // Remove single quotes and whitespace + item.erase(std::remove(item.begin(), item.end(), '\''), item.end()); + item.erase(0, item.find_first_not_of(" \t")); + item.erase(item.find_last_not_of(" \t") + 1); + + if (item.find("Region") != std::string::npos) { + item = "../" + item; + } + + fieldList.append(word(item)); + } + + Info << "Parsed fields:" << nl; + forAll(fieldList, i) + { + Info << " - " << fieldList[i] << nl; + } + + // Loop for field + snapshot sampling + while (true) + { + // Define strings + string snapshotPath; + + // Get snapshot path + std::cout << "Enter snapshot path: "; + if (!std::getline(std::cin, snapshotPath) || snapshotPath.empty()) + { + break; + } + + // Construct snapshot path + std::string fieldTime = snapshotPath; + + forAll(fieldList, fieldI) + { + volScalarField field + ( + IOobject + ( + fieldList[fieldI], + fieldTime, + mesh(), + IOobject::MUST_READ, + IOobject::NO_WRITE + ), + mesh() + ); + + // Initialize variables to store max and min values and their positions + scalar maxVal = -VGREAT; + scalar minVal = VGREAT; + label maxPos = -1; + label minPos = -1; + + // Iterate over the field to find max and min values and their positions + forAll(field, fieldI) + { + if (field[fieldI] > maxVal) + { + maxVal = field[fieldI]; + maxPos = fieldI; + } + if (field[fieldI] < minVal) + { + minVal = field[fieldI]; + minPos = fieldI; + } + } + + // Convert cell labels to positions + const vector& maxPosition = mesh().C()[maxPos]; + const vector& minPosition = mesh().C()[minPos]; + + // Info output + string snapshotID = snapshotPath.substr(snapshotPath.rfind('/', snapshotPath.rfind('/') - 1) + 1); + Info << "Minimum: " << snapshotID << ", location = " << minPosition << ", " << fieldList[fieldI] << " = " << minVal << endl; + Info << "Maximum: " << snapshotID << ", location = " << maxPosition << ", " << fieldList[fieldI] << " = " << maxVal << endl; + } + } + + #}; +} + +// ************************************************************************* // + diff --git a/tools/linAlg4Foam/OpenFoamFunctionObjects/FOinternalFieldMinMax b/tools/linAlg4Foam/OpenFoamFunctionObjects/FOinternalFieldMinMax new file mode 100644 index 0000000..128f411 --- /dev/null +++ b/tools/linAlg4Foam/OpenFoamFunctionObjects/FOinternalFieldMinMax @@ -0,0 +1,130 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: v2406 | +| \\ / A nd | Website: www.openfoam.com | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ + +// This function object calculates the min/max values and their positions +// in the internal field. +// +// Input arguments: +// - path to the snapshot (time folder) +// - field +// + +internalFieldMinMax +{ + type coded; + libs (utilityFunctionObjects); + region region0; + + name internalFieldMinMax; + + codeWrite + #{ + // Read the list of fields in python format (e.g., ['T', 'DT']) + string inputLine; + std::cout << "Enter list of fields: "; + std::getline(std::cin, inputLine); + + // Strip square brackets + if (!inputLine.empty() && inputLine.front() == '[') inputLine.erase(0, 1); + if (!inputLine.empty() && inputLine.back() == ']') inputLine.pop_back(); + + // Split and clean elements + wordList fieldList; + std::stringstream ss(inputLine); + std::string item; + + while (std::getline(ss, item, ',')) { + // Remove single quotes and whitespace + item.erase(std::remove(item.begin(), item.end(), '\''), item.end()); + item.erase(0, item.find_first_not_of(" \t")); + item.erase(item.find_last_not_of(" \t") + 1); + + if (item.find("Region") != std::string::npos) { + item = "../" + item; + } + + fieldList.append(word(item)); + } + + Info << "Parsed fields:" << nl; + forAll(fieldList, i) + { + Info << " - " << fieldList[i] << nl; + } + + // Loop for field + snapshot sampling + while (true) + { + // Define strings + string snapshotPath; + + // Get snapshot path + std::cout << "Enter snapshot path: "; + if (!std::getline(std::cin, snapshotPath) || snapshotPath.empty()) + { + break; + } + + // Construct snapshot path + std::string fieldTime = snapshotPath; + + forAll(fieldList, fieldI) + { + volScalarField field + ( + IOobject + ( + fieldList[fieldI], + fieldTime, + mesh(), + IOobject::MUST_READ, + IOobject::NO_WRITE + ), + mesh() + ); + + // Access only the internal field + const scalarField& internalField = field.internalField(); + + // Initialize variables to store max and min values and their positions + scalar maxVal = -VGREAT; + scalar minVal = VGREAT; + label maxPos = -1; + label minPos = -1; + + // Iterate over the internal field to find max and min values and their positions + forAll(internalField, fieldI) + { + if (internalField[fieldI] > maxVal) + { + maxVal = internalField[fieldI]; + maxPos = fieldI; + } + if (internalField[fieldI] < minVal) + { + minVal = internalField[fieldI]; + minPos = fieldI; + } + } + + // Convert cell labels to positions + const vector& maxPosition = mesh().C()[maxPos]; + const vector& minPosition = mesh().C()[minPos]; + + // Info output + string snapshotID = snapshotPath.substr(snapshotPath.rfind('/', snapshotPath.rfind('/') - 1) + 1); + Info << "Minimum: " << snapshotID << ", location = " << minPosition << ", " << fieldList[fieldI] << " = " << minVal << endl; + Info << "Maximum: " << snapshotID << ", location = " << maxPosition << ", " << fieldList[fieldI] << " = " << maxVal << endl; + } + } + + #}; +} + +// ************************************************************************* // + diff --git a/tools/linAlg4Foam/OpenFoamFunctionObjects/FOlinearCombination b/tools/linAlg4Foam/OpenFoamFunctionObjects/FOlinearCombination new file mode 100644 index 0000000..ca91e9a --- /dev/null +++ b/tools/linAlg4Foam/OpenFoamFunctionObjects/FOlinearCombination @@ -0,0 +1,269 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ______ _ __ ______ | +| / ____/ ___ / | / / / ____/ ____ ____ _ ____ ___ | +| / / __ / _ \ / |/ / ______ / /_ / __ \ / __ `/ / __ `__ \ | +| / /_/ / / __/ / /| / /_____/ / __/ / /_/ // /_/ / / / / / / / | +| \____/ \___/ /_/ |_/ /_/ \____/ \__,_/ /_/ /_/ /_/ | +| Copyright (C) 2015 - 2022 EPFL | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object functions; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +// This function object takes a list of scalar and a list of snapshots paths as +// input arguments, and returns a sum of fields, where each of the field is +// multiplied by appropriate. + +linearCombination +{ + type coded; + libs (utilityFunctionObjects); + region region0; + + name linearCombination; + + codeWrite + #{ + + // Read a list of fields + string fieldsList; + cout << "Enter list of fields: "; + std::getline(std::cin, fieldsList); + + // Remove square brackets if they exist + if (!fieldsList.empty() && fieldsList.front() == '[') fieldsList.erase(0, 1); + if (!fieldsList.empty() && fieldsList.back() == ']') fieldsList.pop_back(); + + // Prepare stream + std::stringstream fieldsStream(fieldsList); + std::string fieldToken; + List fieldNames; + while (std::getline(fieldsStream, fieldToken, ',')) { + // Trim whitespace + fieldToken.erase(0, fieldToken.find_first_not_of(" \t")); + fieldToken.erase(fieldToken.find_last_not_of(" \t") + 1); + // Remove quotes + if (!fieldToken.empty() && (fieldToken.front() == '\'' || fieldToken.front() == '"')) fieldToken.erase(0, 1); + if (!fieldToken.empty() && (fieldToken.back() == '\'' || fieldToken.back() == '"')) fieldToken.pop_back(); + + if (fieldToken.find("Region") != std::string::npos) { + fieldToken = "../" + fieldToken; + } + + fieldNames.append(word(fieldToken)); + } + + // Print the field names + Info << "Parsed field names: " << nl; + forAll(fieldNames, field) { + Info << " " << fieldNames[field] << nl; + } + + while (true) + { + // Read a list of coefficients in a python format separated by a comma (e.g., [0.5, 0.5, 0.5, 0.5]) + string inputCoefList; + std::cout << "Enter list of coefficients: "; + if (!std::getline(std::cin, inputCoefList) || inputCoefList.empty()) + { + break; + } + + // Remove square brackets if they exist + if (!inputCoefList.empty() && inputCoefList.front() == '[') inputCoefList.erase(0, 1); + if (!inputCoefList.empty() && inputCoefList.back() == ']') inputCoefList.pop_back(); + + // Create a list of coefficients + std::stringstream ss(inputCoefList); + std::string token; + List coefficients; + + while (std::getline(ss, token, ',')) { + // Trim whitespace + token.erase(0, token.find_first_not_of(" \t")); + token.erase(token.find_last_not_of(" \t") + 1); + + try { + scalar value = std::stod(token); + coefficients.append(value); + } catch (...) { + Info << "Invalid number: " << token << nl; + } + } + + // Print the coefficients + Info << "Parsed coefficients: " << nl; + forAll(coefficients, i) { + Info << " " << coefficients[i] << nl; + } + + // Read a list of snapshots in a python format (separated by a comma) + string inputSnapshotsList; + std::cout << "Enter list of snapshots: "; + + if (!std::getline(std::cin, inputSnapshotsList) || inputSnapshotsList.empty()) + { + break; + } + + // Remove square brackets if they exist + if (!inputSnapshotsList.empty() && inputSnapshotsList.front() == '[') inputSnapshotsList.erase(0, 1); + if (!inputSnapshotsList.empty() && inputSnapshotsList.back() == ']') inputSnapshotsList.pop_back(); + + // Prepare stream + std::stringstream snapsStream(inputSnapshotsList); + std::string snapsToken; + List snapshotsList; + + while (std::getline(snapsStream, snapsToken, ',')) { + // Trim whitespace + snapsToken.erase(0, snapsToken.find_first_not_of(" \t")); + snapsToken.erase(snapsToken.find_last_not_of(" \t") + 1); + // Remove quotes + if (!snapsToken.empty() && (snapsToken.front() == '\'' || snapsToken.front() == '"')) snapsToken.erase(0, 1); + if (!snapsToken.empty() && (snapsToken.back() == '\'' || snapsToken.back() == '"')) snapsToken.pop_back(); + + snapshotsList.append(word(snapsToken)); + } + + // Print the snapshots + Info << "Parsed snapshots: " << nl; + forAll(snapshotsList, snapshot) { + Info << " " << snapshotsList[snapshot] << nl; + } + + // Check if the number of coefficients matches the number of snapshots + if (coefficients.size() != snapshotsList.size()) + { + FatalErrorInFunction + << "The number of coefficients (" << coefficients.size() << ") does not match the number of snapshots (" << snapshotsList.size() << ")" + << exit(FatalError); + } + + // Read the output location + string outputLocation; + std::cout << "Enter output location: "; + std::getline(std::cin, outputLocation); + + // Remove square brackets + if (!outputLocation.empty() && outputLocation.front() == '[') outputLocation.erase(0, 1); + if (!outputLocation.empty() && outputLocation.back() == ']') outputLocation.pop_back(); + // Remove leading and trailing whitespace + outputLocation.erase(0, outputLocation.find_first_not_of(" \t")); + outputLocation.erase(outputLocation.find_last_not_of(" \t") + 1); + // Remove quotes + if (!outputLocation.empty() && (outputLocation.front() == '\'' || outputLocation.front() == '"')) outputLocation.erase(0, 1); + if (!outputLocation.empty() && (outputLocation.back() == '\'' || outputLocation.back() == '"')) outputLocation.pop_back(); + + // Get the OpenFoam work directory and construct the links + const char* caseDir = "FOAM_CASE"; + const char* caseDirPath = std::getenv(caseDir); + outputLocation = std::string(caseDirPath) + "/" + outputLocation; + + // Make the output directory + if (!isDir(outputLocation)) + { + mkDir(outputLocation); + } + + // Construct path to the first snapshot + std::string fieldTime0 = snapshotsList[0]; + + forAll(fieldNames, field) + { + // Use the first field as a base + volScalarField sumOfFields + ( + IOobject + ( + fieldNames[field], + fieldTime0, + mesh(), + IOobject::MUST_READ, + IOobject::NO_WRITE + ), + mesh() + ); + + // Multiply the field by its coefficient + sumOfFields = coefficients[0] * sumOfFields; + + // Sum the first field with the rest of the fields + if (coefficients.size() > 1) + { + for (label i = 1; i < coefficients.size(); ++i) + { + // Construct path to the snapshot + std::string fieldTime = snapshotsList[i]; + + volScalarField originalField + ( + IOobject + ( + fieldNames[field], + fieldTime, + mesh(), + IOobject::MUST_READ, + IOobject::NO_WRITE + ), + mesh() + ); + + // Multiply the field by its coefficient and add to the first field + originalField = coefficients[i] * originalField; + sumOfFields += originalField; + } + } + /* + // Rename the sum of fields + volScalarField renamedField + ( + IOobject + ( + "linearCombination_" + fieldNames[field], + mesh().time().timeName(), + mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + sumOfFields + ); + + */ + // Extract just the field name without the region path + word cleanFieldName = fileName(fieldNames[field]).name(); + + // Then use cleanFieldName instead of fieldNames[field] + volScalarField renamedField + ( + IOobject + ( + cleanFieldName, // Use cleanFieldName here + mesh().time().timeName(), + mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + sumOfFields + ); + // Result is printed into 'constant' + renamedField.write(); + + // Move the result to the output location + std::string moveCommand = "mv " + renamedField.path() + "/" + cleanFieldName + " " + outputLocation; + system(moveCommand); + + } + + } + + #}; +} + + +// ************************************************************************* // diff --git a/tools/linAlg4Foam/OpenFoamFunctionObjects/FOsample b/tools/linAlg4Foam/OpenFoamFunctionObjects/FOsample new file mode 100644 index 0000000..2d5e648 --- /dev/null +++ b/tools/linAlg4Foam/OpenFoamFunctionObjects/FOsample @@ -0,0 +1,185 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ______ _ __ ______ | +| / ____/ ___ / | / / / ____/ ____ ____ _ ____ ___ | +| / / __ / _ \ / |/ / ______ / /_ / __ \ / __ `/ / __ `__ \ | +| / /_/ / / __/ / /| / /_____/ / __/ / /_/ // /_/ / / / / / / / | +| \____/ \___/ /_/ |_/ /_/ \____/ \__,_/ /_/ /_/ /_/ | +| Copyright (C) 2015 - 2022 EPFL | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object functions; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +// Reads the value of a field in a given point. It is designed to read multiple +// cases, as well as fields and also points. +// +// Input arguments: +// - path to the snapshot (time folder) +// - field +// - point +// + +sample +{ + type coded; + libs (utilityFunctionObjects); + region region0; + + name sample; + + codeWrite + #{ + // Read a list of points in python format (e.g., [(1.0, 1.0, 1.0), (1.5, 1.7, 3.9)]) + std::string inputPoints; + std::cout << "Enter list of points: "; + std::getline(std::cin, inputPoints); + + // Remove surrounding square brackets if present + if (!inputPoints.empty() && inputPoints.front() == '[') inputPoints.erase(0, 1); + if (!inputPoints.empty() && inputPoints.back() == ']') inputPoints.pop_back(); + + // Prepare container for points + List probePoints; + + // Tokenize by ')', each representing the end of a point tuple + std::stringstream pointStream(inputPoints); + std::string segment; + while (std::getline(pointStream, segment, ')')) + { + // Skip empty segments + if (segment.find(',') == std::string::npos) continue; + + // Remove leading '(' or ',' if present + size_t start = segment.find('('); + if (start != std::string::npos) + segment = segment.substr(start + 1); + else + segment.erase(0, segment.find_first_not_of(" ,")); + + // Split segment into coordinates + std::stringstream coordStream(segment); + std::string value; + scalar coords[3]; + int i = 0; + + while (std::getline(coordStream, value, ',') && i < 3) + { + coords[i++] = std::stof(value); + } + + if (i != 3) + { + Info << "Invalid point: (" << segment << ")" << nl; + continue; + } + + probePoints.append(point(coords[0], coords[1], coords[2])); + } + + // Print parsed points + Info << "Parsed points:" << nl; + forAll(probePoints, pointI) + { + Info << " - " << probePoints[pointI] << nl; + } + + // Read the list of fields in python format (e.g., ['T', 'DT']) + std::string inputLine; + std::cout << "Enter list of fields: "; + std::getline(std::cin, inputLine); + + // Strip square brackets + if (!inputLine.empty() && inputLine.front() == '[') inputLine.erase(0, 1); + if (!inputLine.empty() && inputLine.back() == ']') inputLine.pop_back(); + + // Split and clean elements + wordList fieldList; + std::stringstream fieldStream(inputLine); + std::string item; + + while (std::getline(fieldStream, item, ',')) { + // Remove single quotes and whitespace + item.erase(std::remove(item.begin(), item.end(), '\''), item.end()); + item.erase(0, item.find_first_not_of(" \t")); + item.erase(item.find_last_not_of(" \t") + 1); + + // Prepend "../" only if the field contains "Region" + if (item.find("Region") != std::string::npos) { + item = "../" + item; + } + + fieldList.append(word(item)); + } + + Info << "Parsed fields:" << nl; + forAll(fieldList, i) + { + Info << " - " << fieldList[i] << nl; + } + + // Loop for field + snapshot sampling + while (true) + { + // Define strings + string snapshotPath; + + // Get snapshot path + std::cout << "Enter snapshot path: "; + if (!std::getline(std::cin, snapshotPath) || snapshotPath.empty()) + { + break; + } + + // Construct snapshot path + std::string fieldTime = snapshotPath; + + forAll(fieldList, fieldI) + { + volScalarField field + ( + IOobject + ( + fieldList[fieldI], + fieldTime, + mesh(), + IOobject::MUST_READ, + IOobject::NO_WRITE + ), + mesh() + ); + + // Access only the internal field + // const scalarField& internalField = field.internalField(); + + forAll(probePoints, pointI) + { + label cellID = mesh().findCell(probePoints[pointI]); + + // Chech the point is in the domain + if (cellID == -1) + { + Info << "Point " << probePoints[pointI] << " is outside the mesh." << nl; + } + else + { + // Read the field value within that cell + // scalar fieldValue = internalField[cellID]; + scalar fieldValue = field[cellID]; + + // Info output + string snapshotID = snapshotPath.substr(snapshotPath.rfind('/', snapshotPath.rfind('/') - 1) + 1); + Info << snapshotID << ", point = " << probePoints[pointI] << ", " << fieldList[fieldI] << " = " << fieldValue << endl; + } + } + } + } + + #}; +} + +// ************************************************************************* // \ No newline at end of file diff --git a/tools/linAlg4Foam/OpenFoamFunctionObjects/FOsensors b/tools/linAlg4Foam/OpenFoamFunctionObjects/FOsensors new file mode 100644 index 0000000..3439f60 --- /dev/null +++ b/tools/linAlg4Foam/OpenFoamFunctionObjects/FOsensors @@ -0,0 +1,149 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ______ _ __ ______ | +| / ____/ ___ / | / / / ____/ ____ ____ _ ____ ___ | +| / / __ / _ \ / |/ / ______ / /_ / __ \ / __ `/ / __ `__ \ | +| / /_/ / / __/ / /| / /_____/ / __/ / /_/ // /_/ / / / / / / / | +| \____/ \___/ /_/ |_/ /_/ \____/ \__,_/ /_/ /_/ /_/ | +| Copyright (C) 2015 - 2022 EPFL | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object functions; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +// Calculates field parameters over a volume defined in cellSet (topoSet) +// - min - uses gMin in the cellSet +// - max - uses gMin in the cellSet +// - average - uses integral over the cellSet + +sensors +{ + type coded; + libs (utilityFunctionObjects); + region region0; + + name sensor; + + codeWrite + #{ + // Get a list with all the sensors named "sensor" + List sensors; + forAll(mesh().cellZones().names(), zoneI) + { + if (mesh().cellZones().names()[zoneI].find("sensor") != string::npos) + { + sensors.append(mesh().cellZones().names()[zoneI]); + } + } + + // Read the list of fields in python format (e.g., ['T', 'DT']) + std::string inputLine; + std::cout << "Enter list of fields: "; + std::getline(std::cin, inputLine); + + // Strip square brackets + if (!inputLine.empty() && inputLine.front() == '[') inputLine.erase(0, 1); + if (!inputLine.empty() && inputLine.back() == ']') inputLine.pop_back(); + + // Split and clean elements + wordList fieldList; + std::stringstream ss(inputLine); + std::string item; + + while (std::getline(ss, item, ',')) { + // Remove single quotes and whitespace + item.erase(std::remove(item.begin(), item.end(), '\''), item.end()); + item.erase(0, item.find_first_not_of(" \t")); + item.erase(item.find_last_not_of(" \t") + 1); + + // Prepend "../" only if the field contains "Region" + if (item.find("Region") != std::string::npos) { + item = "../" + item; + } + + fieldList.append(word(item)); + } + + Info << "Parsed fields:" << nl; + forAll(fieldList, i) + { + Info << " - " << fieldList[i] << nl; + } + + // Loop for field + snapshot sampling + while (true) + { + // Define strings + string snapshotPath; + + // Get snapshot path + std::cout << "Enter snapshot path: "; + if (!std::getline(std::cin, snapshotPath) || snapshotPath.empty()) + { + break; + } + + // Construct snapshot path + std::string fieldTime = snapshotPath; + + forAll(fieldList, fieldI) + { + volScalarField field + ( + IOobject + ( + fieldList[fieldI], + fieldTime, + mesh(), + IOobject::MUST_READ, + IOobject::NO_WRITE + ), + mesh() + ); + + // Find field values for appropriate regions and calculate needed values + for (label j = 0; j < sensors.size(); ++j) + { + // Find appropriate cell set zone + const labelList& cellIDs = mesh().cellZones()[sensors[j]]; + + // Get cell volumes and field values inside the cell set + double cellSetVolume = 0; + double cellSetDomainIntegral = 0; + + scalarField cellSetField(cellIDs.size()); + + for (label i = 0; i < cellIDs.size(); ++i) + { + label cellID = cellIDs[i]; + cellSetField[i] = field[cellID]; + cellSetVolume = cellSetVolume + mesh().V()[cellID]; + + // Calculate domain integral + cellSetDomainIntegral = cellSetDomainIntegral + mesh().V()[cellID] * field[cellID]; + } + + // Extract min, max and average values + double minValue = gMin(cellSetField); + double maxValue = gMax(cellSetField); + double averageValue = cellSetDomainIntegral/cellSetVolume; + // Info << "Total volume " << cellSetVolume << endl; + // Info << "Domain integral = " << cellSetDomainIntegral << endl; + + // Info output + string snapshotID = snapshotPath.substr(snapshotPath.rfind('/', snapshotPath.rfind('/') - 1) + 1); + Info << snapshotID << ", " << sensors[j] << ", " << fieldList[fieldI] <<", min = " << minValue << ", max = " << maxValue << ", average = " << averageValue << endl; + } + + } + } + + #}; +} + + +// ************************************************************************* // \ No newline at end of file diff --git a/tools/linAlg4Foam/OpenFoamFunctionObjects/controlDict b/tools/linAlg4Foam/OpenFoamFunctionObjects/controlDict new file mode 100644 index 0000000..9a0a100 --- /dev/null +++ b/tools/linAlg4Foam/OpenFoamFunctionObjects/controlDict @@ -0,0 +1,47 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: v2406 | +| \\ / A nd | Website: www.openfoam.com | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object controlDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +application laplacianFoam; + +startFrom startTime; + +startTime 0; + +stopAt endTime; + +endTime 1; + +deltaT 1; + +writeControl runTime; + +writeInterval 5; + +purgeWrite 0; + +writeFormat ascii; + +writePrecision 6; + +writeCompression off; + +timeFormat general; + +timePrecision 6; + +runTimeModifiable true; + +// ************************************************************************* // diff --git a/tools/linAlg4Foam/OpenFoamFunctionObjects/topoSetDict b/tools/linAlg4Foam/OpenFoamFunctionObjects/topoSetDict new file mode 100644 index 0000000..67241f5 --- /dev/null +++ b/tools/linAlg4Foam/OpenFoamFunctionObjects/topoSetDict @@ -0,0 +1,23 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: v2312 | +| \\ / A nd | Website: www.openfoam.com | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object topoSetDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +actions +( + +); + + +// ************************************************************************* // diff --git a/tools/linAlg4Foam/README.md b/tools/linAlg4Foam/README.md new file mode 100644 index 0000000..0fdaa2b --- /dev/null +++ b/tools/linAlg4Foam/README.md @@ -0,0 +1,231 @@ +# Linear algebra for OpenFOAM + +The 'linAlg4Foam.py' contains a collection of field operations which are done inside the +OpenFOAM environment and can be directly called from python. + +## Installation + +Download zip or tar.gz file and extract it. + + +## Usage + +Import the linAlg4Foam into your script as: +``` +from linAlg4Foam import linAlg4Foam +``` + +Then, refer to the utilities as: +``` +linAlg4Foam.utility(arguments) +``` + +If the default directory structure generated by the `snapshot_manager.py` is changed, relative path from OpenFOAM virtual folder to the snapshots location should be adjusted accordingly. The `relative_path` variable is defined in 'basicFunction.py' as: +``` +relative_path = "../symlinked_cases/" +``` + +Check test.py for examples off all implemented functions. + + +## Functions + +### Basic utilities + +#### Extracting minimal and maximal values +`fieldMinMax` and `internalFieldMinMax` return the minimal and maximal values of the fields along with the location. + +Input example +``` +from linAlg4Foam import linAlg4Foam + +# Path to virtual OpenFoam directory +virtual_OF_dir = "OpenFoamTestCase/virtual_OF" + +snapshots_list =[ + '1/100', + '3/5' +] + +fields = ['T', 'DT'] + +## Returning minimal and maximal values of a field +minima, maxima = linAlg4Foam.fieldMinMax(virtual_OF_dir, snapshots_list, fields) +for item in minima: + print(item) +for item in maxima: + print(item) + +## Returning minimal and maximal values of internal field +minima, maxima = linAlg4Foam.internalFieldMinMax(virtual_OF_dir, snapshots_list, fields) +for item in minima: + print(item) +for item in maxima: + print(item) +``` + +Output example +``` +... +{'snapshotID': '1/100', 'field': 'T', 'position': '(0.464286 0.0294118 0.05)', 'value': 300.132} +{'snapshotID': '1/100', 'field': 'DT', 'position': '(0.0357143 0.0294118 0.05)', 'value': 0.000284876} +{'snapshotID': '3/5', 'field': 'T', 'position': '(0.178571 0.0294118 0.05)', 'value': 300.0} +{'snapshotID': '3/5', 'field': 'DT', 'position': '(0.0357143 0.0294118 0.05)', 'value': 0.000143284} +{'snapshotID': '3/200', 'field': 'T', 'position': '(0.0357143 0.5 0.05)', 'value': 270.733} +{'snapshotID': '3/200', 'field': 'DT', 'position': '(0.0357143 0.0294118 0.05)', 'value': 0.000143284} +``` + + +#### Sampling (probing) + +`readFromPoints` reads field values from the set of points, and returns a list of dictionaries. + +Input example +``` +from linAlg4Foam import linAlg4Foam + +# Path to virtual OpenFoam directory +virtual_OF_dir = "OpenFoamTestCase/virtual_OF" + +snapshots_list =[ + '1/100', + '3/5' +] + +points = [ + (0.1, 0.5, 0.05), + (0.2, 0.5, 0.05), + (0.4, 0.5, 0.05) +] + +fields = ['T', 'DT'] + +## Sampling from points +results = linAlg4Foam.readFromPoints(virtual_OF_dir, snapshots_list, points, fields) +for item in results: + print(item) + +``` + + +Output example +``` +... +{'snapshotID': '7/130', 'point': '(0.1 0.5 0.05)', 'field': 'T', 'value': 336.211} +{'snapshotID': '7/130', 'point': '(0.2 0.5 0.05)', 'field': 'T', 'value': 333.071} +{'snapshotID': '7/130', 'point': '(0.4 0.5 0.05)', 'field': 'T', 'value': 307.728} +{'snapshotID': '7/130', 'point': '(0.1 0.5 0.05)', 'field': 'DT', 'value': 0.000275331} +{'snapshotID': '7/130', 'point': '(0.2 0.5 0.05)', 'field': 'DT', 'value': 0.000275331} +{'snapshotID': '7/130', 'point': '(0.4 0.5 0.05)', 'field': 'DT', 'value': 0.000275331} + +``` + +#### Read from a physical sensor + +`physicalSensors` extracts minimal, maximal and average value from a 3D cylindrical part of the domain. This functionality firstly generates cell zones using 'topoSet' and then reads the fields from that specific zones. + +The sensors are listed as a list of dictionaries. Each sensor is defined with two points (on each side of the axis) and a radius. + +Input example +``` +from linAlg4Foam import linAlg4Foam + +# Path to virtual OpenFoam directory +virtual_OF_dir = "OpenFoamTestCase/virtual_OF" + +snapshots_list =[ + '1/100', + '3/5' +] + +fields = ['T', 'DT'] + +sensors = [ + {"point1": (0.0, 0.5, 0.0005), "point2": (0.1, 0.5, 0.0005), "radius": 0.1}, + {"point1": (0.2, 0.5, 0.0005), "point2": (0.3, 0.5, 0.0005), "radius": 0.1} +] + +## Physical sensor +results = linAlg4Foam.physicalSensors(virtual_OF_dir, snapshots_list, sensors, fields) +for item in results: + print(item) + +``` + +Output example +``` +... +{'snapshotID': '7/150', 'sensor': 'sensor0', 'field': 'DT', 'min': 0.000275331, 'max': 0.000275331, 'average': 0.000275331} +{'snapshotID': '7/150', 'sensor': 'sensor1', 'field': 'DT', 'min': 0.000275331, 'max': 0.000275331, 'average': 0.000275331} +{'snapshotID': '7/130', 'sensor': 'sensor0', 'field': 'T', 'min': 319.676, 'max': 320.017, 'average': 319.79} +{'snapshotID': '7/130', 'sensor': 'sensor1', 'field': 'T', 'min': 323.48, 'max': 323.886, 'average': 323.615} +{'snapshotID': '7/130', 'sensor': 'sensor0', 'field': 'DT', 'min': 0.000275331, 'max': 0.000275331, 'average': 0.000275331} +{'snapshotID': '7/130', 'sensor': 'sensor1', 'field': 'DT', 'min': 0.000275331, 'max': 0.000275331, 'average': 0.000275331} + +``` + +### Field Operations + +#### Linear combinations of fields + +`linearCombination` calculates a linear combination of fields as: +$$ +\sum_{i = 1}^N c_i \psi_i +$$ + +The inputs are: +- list of field names (fields), +- list of lists of snapshot paths (list_of_lists_of_snapshots), +- list of lists of coefficients (list_of_lists_of_coefficients), +- list of output directories (output_dirs). + +The output directories should be specified in the same order as the +snapshots, since the function will create a new directory for each +output directory and write new fields there. + + +Input example +``` +from linAlg4Foam import linAlg4Foam + +# Path to virtual OpenFoam directory +virtual_OF_dir = "OpenFoamTestCase/virtual_OF" + +fields = ['T', 'DT'] + +list_of_lists_of_snapshots = [ + ['4/100', '5/100', '6/100'], + ['1/200', '2/200', '3/200'] + ] + +list_of_lists_of_coefficients = [ + [0.5, -0.5, 0.5], + [10.0, -1.0e-4, 0.1] + ] + +output_dirs = ['temp_outputs/0', 'temp_outputs/1'] + +## Linear combination of fields +linAlg4Foam.linearCombination(virtual_OF_dir, list_of_lists_of_snapshots, fields, list_of_lists_of_coefficients, output_dirs) + +``` + + +## Adding new functions + +1) Create new function object ('FOnewFunction') and place it into `OpenFoamFunctionObjects` directory. +2) In linAlg4Foam.py, create a function 'newFunction' and place it under 'class linAlg4Foam'. + +An example of a new function definition (a set of input arguments may vary): +``` +@staticmethod +def newFunction(virtual_OF_dir, snapshots_list, fields, points, list_of_lists_of_coefficients, output_dirs, sensors): + function = 'newFunction' + initialize(virtual_OF_dir, function) + setSensors(virtual_OF_dir, sensors) + execute(virtual_OF_dir=virtual_OF_dir, function=function, snapshots_list=snapshots_list, fields=fields, points=points, list_of_lists_of_coefficients = list_of_lists_of_coefficients, output_dirs = output_dirs) + + return +``` + +3) If an extra prompt is needed, add it into prompt_handlers list in basicFunctions.py. diff --git a/tools/linAlg4Foam/basicFunctions.py b/tools/linAlg4Foam/basicFunctions.py new file mode 100644 index 0000000..3444afe --- /dev/null +++ b/tools/linAlg4Foam/basicFunctions.py @@ -0,0 +1,271 @@ +""" +=============================================================================== + Script Name: basicFunctions.py + Description: This script contains basic file operators which are + used in 'linAlg4Foam.py' + + Author: Rok Krpan (rok.krpan@tamu.edu) + Date: 2025-05-23 + Version: 1.2 + Usage: To properly run this script, one should have OpenFOAM + intalled on the system and OpenFOAM evironment variables + should be properly sourced inside '~/.bashrc' + The functions definitions from this script should + be imported in 'linAlg4Foam.py' as: + 'from basicFunctions import *' + +=============================================================================== +""" + +import gzip +import os +import pexpect +import re +import shutil +import subprocess + +from attr import field +from itertools import combinations + +# Location of OpenFoam function objects +function_objects_dir = os.path.join(os.path.dirname(__file__), 'OpenFoamFunctionObjects') + +# Relative path from OpenFoam virtual folder to the snapshots location +relative_path = "../symlinked_cases/" + +class ExecuteFunctionObject: + def __init__(self, virtual_OF_dir, function, snapshots_list, fields, points = None, list_of_lists_of_coefficients = None, output_dirs = None, sensors = None): + self.virtual_OF_dir = virtual_OF_dir + self.function = function + self.snapshots_list = snapshots_list + self.fields = fields + self.points = points + self.list_of_lists_of_coefficients = list_of_lists_of_coefficients + self.output_dirs = output_dirs + self.sensors = sensors + + self.regions = self._check_regions_() + self.initialize() + if self.function == "sensors": + self.setSensors() + self.execute() + + @staticmethod + def __list_into_string__(value): + return str(value if isinstance(value, list) else [value]) + + @staticmethod + def __open_mesh_file__(mesh_file_path): + gz_path = mesh_file_path + ".gz" + if os.path.isfile(gz_path): + return gzip.open(gz_path, 'rt') + elif os.path.isfile(mesh_file_path): + return open(mesh_file_path, 'r') + else: + raise FileNotFoundError(f"Neither {mesh_file_path} nor {gz_path} found.") + + + # Check the existance of the regions in OpenFOAM case and compare their meshes + def _check_regions_(self): + + if isinstance(self.snapshots_list[0], list): + snapshots_list = [sublist[0] for sublist in self.snapshots_list if sublist] + else: + snapshots_list = self.snapshots_list + + + snapshot_path = os.path.join(self.virtual_OF_dir, relative_path, snapshots_list[0]) + with os.scandir(snapshot_path) as entries: + regions_check = any(entry.is_dir() and "Region" in entry.name for entry in entries) + + # Extract unique regions from fields list + regions = {field.split('/')[0] for field in self.fields if '/' in field} + regions = list(regions) + + + # Check if both cases or regions use the same mesh + if regions_check and len(regions) > 1: + # print(f"Several regions exist: {regions}") + constant_dir = os.path.abspath(os.path.join(snapshot_path, "../constant")) + mesh_files = ['points','faces'] + + for region1, region2 in combinations(regions, 2): + for mesh_file in mesh_files: + mesh_file1 = os.path.join(constant_dir, region1, "polyMesh", mesh_file) + mesh_file2 = os.path.join(constant_dir, region2, "polyMesh", mesh_file) + # print(f"Comparing {mesh_file1} vs {mesh_file2}") + + with self.__open_mesh_file__(mesh_file1) as f1, self.__open_mesh_file__(mesh_file2) as f2: + for i in range(50): # First 50 lines + line1 = f1.readline() + line2 = f2.readline() + + if not line1 and not line2: + break + if not line1 or not line2: + raise ValueError(f"{mesh_file}: {region1} and {region2} differ in file length.") + if line1.strip() != line2.strip(): + raise ValueError(f"{mesh_file}: {region1} and {region2} differ at line {i+1}") + + return regions + + def initialize(self): + os.makedirs(self.virtual_OF_dir + '/system/FOs', exist_ok=True) + os.makedirs(self.virtual_OF_dir + '/logs', exist_ok=True) + + # Copy 'controlDict' from global FO repo + shutil.copy2(function_objects_dir + '/controlDict', self.virtual_OF_dir + '/system/controlDict') + + # Modify 'controlDict' to use the specific function object + command = ["foamDictionary", self.virtual_OF_dir + "/system/controlDict", "-disableFunctionEntries", "-entry", "functions", "-set", "{}" ] + subprocess.run(command, check=True, stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL) + + command = ["foamDictionary", self.virtual_OF_dir + "/system/controlDict", "-disableFunctionEntries", "-entry", "functions/#include", "-set", f"\"FOs/FO{self.function}\"" ] + subprocess.run(command, check=True, stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL) + + # Copy function object if not already there + FoamFile = f"FO{self.function}" + FoamFilePath = os.path.join(self.virtual_OF_dir, 'system', 'FOs') + shutil.copy2(function_objects_dir + f"/FO{self.function}", FoamFilePath) + + # Since the meshes should be the same, there is no difference in the selection of the mesh region. + # If regions exist, select the appropriate region in the function object definition file + if self.regions: + with open(os.path.join(FoamFilePath, FoamFile), "r") as f: + lines = f.readlines() + + with open(os.path.join(FoamFilePath, FoamFile), "w") as f: + for line in lines: + if "region0" in line: + f.write(f" region {self.regions[0]};\n") + else: + f.write(line) + + return + + def setSensors(self): + """ + This function sets cylindrical cell zone sets in 'topoSetDict' to be used as a sensor volumes. + """ + + # Read topoSetDict template + with open(function_objects_dir +"/topoSetDict", "r") as template: + lines = template.readlines() + template_header = lines[:18] + template_footer = lines[19:] + + # Create the topoSetDict in the virtual_OF_dir folder + if self.regions: + output_file = os.path.join(self.virtual_OF_dir, 'system', self.regions[0], 'topoSetDict') + else: + output_file = self.virtual_OF_dir + '/system/topoSetDict' + + # Write topoSetDict + with open(output_file, "w") as f: + f.writelines(template_header) + for i, sensor in enumerate(self.sensors): + # Construct strings for the points (space-separated values inside parentheses) + point1_str = "(" + " ".join(str(x) for x in sensor["point1"]) + ")" + point2_str = "(" + " ".join(str(x) for x in sensor["point2"]) + ")" + + # Construct the dictionary block + block = ( + " {\n" + f" name sensor{i};\n" + " type cellZoneSet;\n" + " action new;\n" + " source cylinderToCell;\n" + f" point1 {point1_str};\n" + f" point2 {point2_str};\n" + f" radius {sensor['radius']};\n" + " }\n" + ) + f.write(block) + f.writelines(template_footer) + + # Execute 'topoSet' utility + # command = ["topoSet", "-case", self.virtual_OF_dir] + command = ["topoSet", "-noFunctionObjects", "-time", "0", "-case", self.virtual_OF_dir] + if self.regions: + command += ["-region", self.regions[0]] + with open(self.virtual_OF_dir + '/logs/topoSet.log', "w") as log_file: + subprocess.run(command, check=True, stdout=log_file, stderr=log_file) + return + + def execute(self): + ## Construct path to snapshots + if isinstance(self.snapshots_list, list) and all(isinstance(sublist, list) for sublist in self.snapshots_list): + snapshots_list = [ + [relative_path + snapshot for snapshot in sublist] + for sublist in self.snapshots_list + ] + else: + snapshots_list = [relative_path + snapshot for snapshot in self.snapshots_list] + + # Construct the command + command = ["postProcess -time '0' -case " + self.virtual_OF_dir] + if self.regions: + command += [" -region ", self.regions[0]] + + # Set log file. + log_path = self.virtual_OF_dir + f"/logs/{self.function}.log" + + # List of prompts and responses + snapshots_iter = iter(snapshots_list) + fields_iter = iter(self.fields) + coefficients_iter = iter(self.list_of_lists_of_coefficients) if self.list_of_lists_of_coefficients is not None else None + output_dirs_iter = iter(self.output_dirs) if self.output_dirs is not None else None + + prompt_handlers = { + "Enter field name:": lambda: next(fields_iter), + "Enter list of fields:": lambda: str(self.fields), + "Enter list of points:": lambda: str(self.points), + "Enter snapshot path: ": lambda: str(next(snapshots_iter)), + "Enter list of snapshots:": lambda: self.__list_into_string__(next(snapshots_iter)), + "Enter list of coefficients:": lambda: self.__list_into_string__(next(coefficients_iter)), + "Enter output location:": lambda: next(output_dirs_iter) + } + + # Spawn the process. + child = pexpect.spawn(" ".join(command), encoding="utf-8", timeout=100) + + with open(log_path, "w") as log_file: + child.logfile = log_file + + try: + while True: + # Build dynamic regex pattern to match any known prompt + prompt_regex = "|".join(map(re.escape, prompt_handlers.keys())) + index = child.expect([prompt_regex, pexpect.EOF]) + + if index == 1: # EOF + break + + matched_prompt = child.match.group(0) + response_func = prompt_handlers.get(matched_prompt) + + if response_func: + try: + response = response_func() + child.sendline(response) + except StopIteration: + # End input if no more data + child.sendeof() + break + + except Exception as e: + print(f"An error occurred: {e}") + + child.close() + + return + + + + + +############################################################################### + +if __name__ == "__main__": + print("Just standing here, doing nothing...") + diff --git a/tools/linAlg4Foam/linAlg4Foam.py b/tools/linAlg4Foam/linAlg4Foam.py new file mode 100644 index 0000000..67d1ef0 --- /dev/null +++ b/tools/linAlg4Foam/linAlg4Foam.py @@ -0,0 +1,324 @@ +""" +=============================================================================== + Script Name: linAlg4Foam.py + Description: This script contains a set of utilites which can be used + in python and perform field operations using OpenFOAM + utilities. It actually acts as an interface which performs + field operation on existing OpenFOAM results using OpenFOAM + function objects. + + Author: Rok Krpan (rok.krpan@tamu.edu) + Date: 2025-05-23 + Version: 1.2 + Usage: To properly run this script, one should have OpenFOAM + installed on the system and OpenFOAM evironment variables + should be properly sourced inside '~/.bashrc'. The functions + defined in this script should be imported as: + 'from linAlg4Foam import linAlg4Foam' + + 'virtual_OF_dir' should be a directiory with "constant" and "system" + folders (with all included files). Namely, the function objects + are executed in this directory and require the basi OpenFOAM + case infomation (polyMesh, fvSystem, fvSolution). + +=============================================================================== +""" + +from .basicFunctions import * + +class linAlg4Foam: + +############################################################################### +### +### Basic utilities +### +############################################################################### + + + @staticmethod + def fieldMinMax(virtual_OF_dir, snapshots_list, fields): + """ + Returns the minimal and maximal values of the fields along with the location. + + Returns: + - output file 'fieldMinMax.log' in virtual_OF_dir/logs/ + - a list of dictionaries + """ + + function = 'fieldMinMax' + ExecuteFunctionObject(virtual_OF_dir=virtual_OF_dir, function=function, snapshots_list=snapshots_list, fields=fields) + + # Extract data from log file + log_path = virtual_OF_dir + f"/logs/{function}.log" + minima = [] + maxima = [] + minimum_pattern = r'Minimum: "([^"]+)", location = \(([^)]+)\), ([\w/\.]+) = ([\d.+-eE]+)' + maximum_pattern = r'Maximum: "([^"]+)", location = \(([^)]+)\), ([\w/\.]+) = ([\d.+-eE]+)' + + with open(log_path, "r") as file: + for line in file: + minimum_match = re.search(minimum_pattern, line) + if minimum_match: + min_snapshotID = minimum_match.group(1) + min_location = f"({minimum_match.group(2)})" + min_field = minimum_match.group(3) + min_value = float(minimum_match.group(4)) + + minima.append({ + "snapshotID": min_snapshotID, + "field": min_field, + "position": min_location, + "value": min_value + }) + + maximum_match = re.search(maximum_pattern, line) + if maximum_match: + max_snapshotID = maximum_match.group(1) + max_location = f"({maximum_match.group(2)})" + max_field = maximum_match.group(3) + max_value = float(maximum_match.group(4)) + + maxima.append({ + "snapshotID": max_snapshotID, + "field": max_field, + "position": max_location, + "value": max_value + }) + + return minima,maxima + + + @staticmethod + def internalFieldMinMax(virtual_OF_dir, snapshots_list, fields): + """ + Returns the minimal and maximal values of the internal fields along + with the location. + + Returns: + - output file 'internalFieldMinMax.log' in virtual_OF_dir/logs/ + - a list of dictionaries + """ + + function = 'internalFieldMinMax' + ExecuteFunctionObject(virtual_OF_dir=virtual_OF_dir, function=function, snapshots_list=snapshots_list, fields=fields) + + # Extract data from log file + log_path = os.path.join(virtual_OF_dir, 'logs', f"{function}.log") + # double dictionary of list initialization + minima = {} + maxima = {} + for snap in snapshots_list: + minima[snap] = {} + maxima[snap] = {} + for field in fields: + minima[snap][field] = [] + maxima[snap][field] = [] + + + minimum_pattern = r'Minimum: "([^"]+)", location = \(([^)]+)\), ([\w/\.]+) = ([\d.+-eE]+)' + maximum_pattern = r'Maximum: "([^"]+)", location = \(([^)]+)\), ([\w/\.]+) = ([\d.+-eE]+)' + + with open(log_path, "r") as file: + for line in file: + minimum_match = re.search(minimum_pattern, line) + if minimum_match: + min_snapshotID = minimum_match.group(1) + min_location = f"({minimum_match.group(2)})" + min_field = minimum_match.group(3) + min_value = float(minimum_match.group(4)) + field_name = min_field.removeprefix('../') + + minima[min_snapshotID][field_name].extend([min_location, min_value]) + + maximum_match = re.search(maximum_pattern, line) + if maximum_match: + max_snapshotID = maximum_match.group(1) + max_location = f"({maximum_match.group(2)})" + max_field = maximum_match.group(3) + max_value = float(maximum_match.group(4)) + field_name = max_field.removeprefix('../') + + maxima[max_snapshotID][field_name].extend([max_location, max_value]) + + return minima,maxima + + @staticmethod + def L2norm(virtual_OF_dir, snapshots_list, fields): + """ + Returns the L2norm of the field. + + Returns: + - output file 'L2norm.log' in virtual_OF_dir/logs/ + - a list of dictionaries with L2norm values + """ + + function = 'L2norm' + ExecuteFunctionObject(virtual_OF_dir=virtual_OF_dir, function=function, snapshots_list=snapshots_list, fields=fields) + + # Extract data from log file + log_path = os.path.join(virtual_OF_dir, 'logs', f"{function}.log") + + # double dictionary of list initialization + L2norm = {} + for snap in snapshots_list: + L2norm[snap] = {} + for field in fields: + L2norm[snap][field] = [] + pattern = r'L2norm: "([^"]+)", (\.\./)?([\w/]+) = ([\d.eE+-]+)' + + with open(log_path, "r") as file: + for line in file: + match = re.search(pattern, line) + if match: + snapshotID = match.group(1) + field = match.group(3) + value = float(match.group(4)) + L2norm[snapshotID][field].extend([value]) + + return L2norm + + @staticmethod + def readFromPoints(virtual_OF_dir, snapshots_list, fields, points): + """ + This function reads field values from the set of points, which are specified as: + points = [ + (0.1, 0.5, 0.05), + (0.2, 0.5, 0.05), + (0.4, 0.5, 0.05) + ] + + Returns: + - output file 'sample.log' in virtual_OF_dir/logs/ + - a list of dictionaries + """ + + function = 'sample' + ExecuteFunctionObject(virtual_OF_dir=virtual_OF_dir, function=function, snapshots_list=snapshots_list, fields=fields, points=points) + + # Extract data from log file + log_path = virtual_OF_dir + f"/logs/{function}.log" + + # double dictionary of list initialization + dirac_delta_evaluation = {} + for snap in snapshots_list: + dirac_delta_evaluation[snap] = {} + for field in fields: + dirac_delta_evaluation[snap][field] = [] + + results = [] + pattern = r'"([^"]+)", point = \(([^)]+)\), ([\w/\.]+) = ([\d.+-eE]+)' + with open(log_path, "r") as file: + for line in file: + match = re.search(pattern, line) + if match: + snapshotID = match.group(1) + point = f"({match.group(2)})" + field = match.group(3).lstrip("../") + value = float(match.group(4)) + dirac_delta_evaluation[snapshotID][field].extend([point, value]) + + return dirac_delta_evaluation + + + + @staticmethod + def physicalSensors(virtual_OF_dir, snapshots_list, fields, sensors): + """ + Extracts minimal, maximal and average value from a 3D cylindrical + part of the domain. This functionality firstly generates cell zones + using 'topoSet' and then reads the fields from that specific zones. + + The sensors are listed as a list of dictionaries. Each sensor is + defined with two points (on each side of the axis) and a radius: + + sensors = [ + {"point1": (0.0, 0.5, 0.0005), "point2": (0.1, 0.5, 0.0005), "radius": 0.1}, + {"point1": (0.2, 0.5, 0.0005), "point2": (0.3, 0.5, 0.0005), "radius": 0.1} + ] + + Returns: + - output file 'physicalSensors.log' in virtual_OF_dir/logs/ + - a list of dictionaries + """ + + function = 'sensors' + ExecuteFunctionObject(virtual_OF_dir=virtual_OF_dir, function=function, snapshots_list=snapshots_list, fields=fields, sensors=sensors) + + # Extract the data from log + log_path = virtual_OF_dir + f"/logs/{function}.log" + results = [] + pattern = r'"([^"]+)",\s*(sensor\d+),\s*([\w/\.]+),\s*min = ([\d.+-eE]+),\s*max = ([\d.+-eE]+),\s*average = ([\d.+-eE]+)' + + with open(log_path, "r") as file: + for line in file: + match = re.search(pattern, line) + if match: + snapshotID = match.group(1) + sensorID = match.group(2) + field = match.group(3).lstrip("../") + min_val = float(match.group(4)) + max_val = float(match.group(5)) + avg_val = float(match.group(6)) + + results.append({ + "snapshotID": snapshotID, + "sensor": sensorID, + "field": field, + "min": min_val, + "max": max_val, + "average": avg_val + }) + + return results + + +############################################################################### +### +### Linear field algebra +### +############################################################################### + + @staticmethod + def linearCombination(virtual_OF_dir, list_of_lists_of_snapshots, fields, list_of_lists_of_coefficients, output_dirs): + """ + Calculates a linear combination of fields: + $$ + \sum_{i = 1}^N c_i \psi_i + $$ + where $c_i$ are the coefficients and $\psi_i$ are the fields. + + The inputs are: + - list of field names (fields), + - list of lists of snapshot paths (snapshots_list), + - list of lists of coefficients (coefficients), + - list of output directories (output_dirs). + + The output directories should be specified in the same order as the + snapshots, since the function will create a new directory for each + output directory and write new fields there. + + Returns: + - output file 'linearCombination.log' in virtual_OF_dir/logs/ + - new fields in output_dirs + """ + + # Check that both lists are lists of lists + if not (isinstance(list_of_lists_of_snapshots, list) and all(isinstance(elem, list) for elem in list_of_lists_of_snapshots)): + print("Error: 'snapshots_list' it is not a list of lists.") + exit(1) + if not (isinstance(list_of_lists_of_coefficients, list) and all(isinstance(elem, list) for elem in list_of_lists_of_coefficients)): + print("Error: 'coefficients' it is not a list of lists.") + exit(1) + + # Check that the lengths of the lists of lists are the same + assert len(list_of_lists_of_snapshots) == len(list_of_lists_of_coefficients) and all(len(s) == len(c) for s, c in zip(list_of_lists_of_snapshots, list_of_lists_of_coefficients)), "Lists are inconsistent" + assert len(list_of_lists_of_snapshots) == len(output_dirs), "Number of output directories must match the number of snapshot lists!" + + function = 'linearCombination' + ExecuteFunctionObject(virtual_OF_dir=virtual_OF_dir, function=function, snapshots_list=list_of_lists_of_snapshots, fields=fields, list_of_lists_of_coefficients=list_of_lists_of_coefficients, output_dirs=output_dirs) + + +############################################################################### + +if __name__ == "__main__": + print("Just standing here, doing nothing...") diff --git a/tools/linAlg4Foam/test.py b/tools/linAlg4Foam/test.py new file mode 100644 index 0000000..348e496 --- /dev/null +++ b/tools/linAlg4Foam/test.py @@ -0,0 +1,86 @@ +from linAlg4Foam import linAlg4Foam + +# Path to virtual OpenFoam directory +virtual_OF_dir = "OpenFoamTestCase/laplacianFoam/virtual_OF" + +#To get value of a scalar/vector field at a specific point inside the domain we need to pass the points to the functions. This is an example how the points can be formed. +# the first, second, third elements of a tuple respectively notifies x, y, and z coordinates of a point +points = [ + (0.1, 0.5, 0.05), + (0.2, 0.5, 0.05), + (0.4, 0.5, 0.05) +] + +#We specify the fields for which we want to do operations +fields = ['T', 'DT'] + +#The portion of string before "/" is the case name. The portion of string after "/" is time name. List of lists of snapshots is for linear combinations (linAlg4Foam.linearCombination). For each list we have to have a corresponding +# list of coefficients and a a output directory. +list_of_lists_of_snapshots = [ + ['4/100', '5/100', '6/100'], + ['1/200', '2/200', '3/200'] +] + +# This coefficients are random. We have put it here to show how the list of lists of coefficients for linearn combination +# can be passed to the linAlg4Foam.linearCombination function. +list_of_lists_of_coefficients = [[0.5, -0.5, 0.5], [10.0, -1.0e-4, 0.1]] +output_dirs = ['temp_outputs/0', 'temp_outputs/1'] + +#Snapshot list for different operations +snapshots_list =[ +'3/5', '3/10', '3/15', '3/20', '3/25', '3/30', '3/35', '3/40', '3/45', '3/50' +] + +#This variable is defined for the linAlg4Foam.physicalSensors function. We define two points and a radius for a physical sensor. +sensors = [ + {"point1": (0.0, 0.5, 0.0005), "point2": (0.1, 0.5, 0.0005), "radius": 0.1}, + {"point1": (0.2, 0.5, 0.0005), "point2": (0.3, 0.5, 0.0005), "radius": 0.1} +] + +## Evaluating fields at the points +print("\Evaluating fields at the points...") +results = linAlg4Foam.readFromPoints(virtual_OF_dir, snapshots_list, fields, points) +#this shows how the results can be shown. +print("Results from the evaluation of linAlg4Foam.readFromPoints") +for item in results: + print(item) + +## Evaluating fields at the physical sensors +print("\nEvaluating at the physical sensors...") +results = linAlg4Foam.physicalSensors(virtual_OF_dir, snapshots_list, fields, sensors) +print("Results from the evaluation of linAlg4Foam.physicalSensors") +for item in results: + print(item) + +## Returning minimal and maximal values of a field +print("\nSearching field min and max...") +minima, maxima = linAlg4Foam.fieldMinMax(virtual_OF_dir, snapshots_list, fields) + +print("Results from the evaluation of linAlg4Foam.fieldMinMax") + +print("minima:\n") +for item in minima: + print(item) + +print("maxima:\n") +for item in maxima: + print(item) + +## Returning minimal and maximal values of internal field +# internal field means fields inside the domain (excluding the boundary) +# linAlg4Foam.internalFieldMinMax determines field min-max inside the domain +print("\nSearching internal field min and max...") +minima, maxima = linAlg4Foam.internalFieldMinMax(virtual_OF_dir, snapshots_list, fields) +print("minima:\n") +for item in minima: + print(item) + +print("maxima:\n") +for item in maxima: + print(item) + +## Linear combination of fields. For each list of fields we get a single linear combination. +print("\nCalculating linear combination of fields...") +linAlg4Foam.linearCombination(virtual_OF_dir, list_of_lists_of_snapshots, fields, list_of_lists_of_coefficients, output_dirs) + + diff --git a/tools/linAlg4Foam/test_multiregion.py b/tools/linAlg4Foam/test_multiregion.py new file mode 100644 index 0000000..4ac9013 --- /dev/null +++ b/tools/linAlg4Foam/test_multiregion.py @@ -0,0 +1,95 @@ +from linAlg4Foam import linAlg4Foam + +# Path to virtual OpenFoam directory +virtual_OF_dir = "OpenFoamTestCase/GeN-Foam/virtual_OF" + +#To get value of a scalar/vector field at a specific point inside the domain we need to pass the points to the functions. This is an example how the points can be formed. +# the first, second, third elements of a tuple respectively notifies x, y, and z coordinates of a point +points = [ + (0, 0.5, -0.6), + (0, 0.5, 0.0), + (0, 0.5, 0.6) +] + +#We specify the fields for which we want to do operations. Notice that we have combined the region with the actual field name to get the constructed field name. +fields = ['fluidRegion/T', 'fluidRegion/p_rgh', 'neutroRegion/flux0'] + +#this field list is for linAlg4Foam.linearCombination +fields_combination = ['fluidRegion/T'] +#The portion of string before "/" is the case name. The portion of string after "/" is time name. List of lists of snapshots is for linear combinations (linAlg4Foam.linearCombination). For each list we have to have a corresponding +# list of coefficients and a a output directory. +list_of_lists_of_snapshots = [ + ['1/1', '1/2', '1/3'], + ['1/4', '1/5'] +] + +# This coefficients are random. We have put it here to show how the list of lists of coefficients for linearn combination +# can be passed to the linAlg4Foam.linearCombination function. +list_of_lists_of_coefficients = [[0.5, -0.5, 0.5], [10.0, -1.0e-4]] +output_dirs = ['temp_outputs/0', 'temp_outputs/1'] + +#Snapshot list for different operations +snapshots_list =[ +'1/1', '1/2', '1/3', '1/4', '1/5' +] + +#This variable is defined for the linAlg4Foam.physicalSensors function. We define two points and a radius for a physical sensor. +sensors = [ + {"point1": (0.0, 0.0, -0.1), "point2": (0.0, 1.9, -0.1), "radius": 0.025}, + {"point1": (0.0, 0.1, 0.1), "point2": (0.0, 2.0, 0.1), "radius": 0.025} +] + +## Evaluating fields at the points +print("\Evaluating fields at the points...") +results = linAlg4Foam.readFromPoints(virtual_OF_dir, snapshots_list, fields, points) +print("Results from the evaluation of linAlg4Foam.readFromPoints.") +for item in results: + print(item) + +## Evaluating fields at the physical sensors +print("\nEvaluating at the physical sensors...") +results = linAlg4Foam.physicalSensors(virtual_OF_dir, snapshots_list, fields, sensors) +print("Results from the evaluation of linAlg4Foam.physicalSensors") +for item in results: + print(item) + +## Returning minimal and maximal values of a field +print("\nSearching field min and max...") +minima, maxima = linAlg4Foam.fieldMinMax(virtual_OF_dir, snapshots_list, fields) +print("Results from the evaluation of linAlg4Foam.fieldMinMax") + +print("minima:\n") +for item in minima: + print(item) + +print("maxima:\n") +for item in maxima: + print(item) + +## Returning minimal and maximal values of internal field +# internal field means fields inside the domain (excluding the boundary) +# linAlg4Foam.internalFieldMinMax determines field min-max inside the domain +print("\nSearching internal field min and max...") +minima, maxima = linAlg4Foam.internalFieldMinMax(virtual_OF_dir, snapshots_list, fields) + +print("Results from the evaluation of linAlg4Foam.internalFieldMinMax") +print("minima:\n") +for item in minima: + print(item) + +print("maxima:\n") +for item in maxima: + print(item) + + +## Returning L2norm of a field +print("\nCalculating L2norm...") +L2norm = linAlg4Foam.L2norm(virtual_OF_dir, snapshots_list, fields) +print("Results from linAlg4Foam.L2norm") +for item in L2norm: + print(item) + +## Linear combination of fields. For each list of fields we get a single linear combination. +print("\nCalculating linear combination of fields...") +linAlg4Foam.linearCombination(virtual_OF_dir, list_of_lists_of_snapshots, fields_combination, list_of_lists_of_coefficients, output_dirs) +