diff --git a/analysis/rev_analysis.cpp b/analysis/rev_analysis.cpp new file mode 100644 index 00000000..c84d0bb2 --- /dev/null +++ b/analysis/rev_analysis.cpp @@ -0,0 +1,864 @@ +/* + Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University + Copyright Equnior ASA + + This file is part of the Open Porous Media project (OPM). + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ +/* + * Representative elementary volume analysis + */ +#include "analysis/rev_analysis.h" +#include "analysis/distance.h" +#include "common/ReadMicroCT.h" +#include "IO/Reader.h" +#include +#include +#include + +ScaLBL_REV::ScaLBL_REV(int RANK, int NP, const Utilities::MPI &COMM) + : average_pore_size(0), rev_x_poro(-1), rev_x_perm(-1), rev_x_tort(-1), + rev_x_surf(-1), size_step(0.06), cc_interval(10), stat_samples(50), gamma(0.1), + tau(1.0), Fx(0), Fy(0), Fz(1.0e-5), + Nx(0), Ny(0), Nz(0), N(0), Np(0), rank(RANK), nprocx(0), nprocy(0), + nprocz(0), nprocs(NP), Lx(0), Ly(0), Lz(0), comm(COMM) {} + +ScaLBL_REV::~ScaLBL_REV() = default; + +void ScaLBL_REV::ReadParams(string filename, bool hydraulic_default) { + + db = std::make_shared(filename); + run_rev = false; + select_subdomains = false; + subdomain_size = {{0, 0, 0}}; + stat_poro_rev = StatisticalREVResult{}; + stat_perm_rev = StatisticalREVResult{}; + if(db->keyExists("REV")){ + run_rev = true; + } + else{ + return; + } + auto mrt_db = db->keyExists("MRT") ? db->getDatabase("MRT") : nullptr; + + size_step = 0.1; + cc_interval = 10; + gamma = 0.1; + stat_samples = 50; + geometrical = true; + hydraulic = hydraulic_default; + standalone_ready = false; + tau = 1.0; + Fx = Fy = 0.0; + Fz = 1.0e-5; + + auto rev_db = db->getDatabase("REV"); + geometrical = rev_db->getWithDefault("Geometrical", true); + hydraulic = rev_db->keyExists("Hydraulic") + ? rev_db->getScalar("Hydraulic") + : hydraulic_default; + seed_anchor = rev_db->getWithDefault("Seed_Anchor", true); + select_subdomains = rev_db->getWithDefault("SubdomainSelect", false); + save_subdomain_raw = rev_db->getWithDefault("save_subdomain_raw", false); + if (rev_db->keyExists("subdomain_size")) { + const std::string message = "REV.subdomain_size requires one or three positive integer voxel dimensions"; + if (!rev_db->isType("subdomain_size")) + throw std::runtime_error(message); + const auto size = rev_db->getVector("subdomain_size"); + if (size.size() != 1 && size.size() != 3) + throw std::runtime_error(message); + for (int d = 0; d < 3; ++d) { + const double value = size[size.size() == 1 ? 0 : d]; + if (!std::isfinite(value) || value < 1.0 || + value > std::numeric_limits::max() || std::floor(value) != value) + throw std::runtime_error(message); + subdomain_size[d] = static_cast(value); + } + } + INSIST(!select_subdomains || comm.getSize() == 1, + "REV.SubdomainSelect requires one MPI process"); + if (rev_db->keyExists("SizeStep")) + size_step = rev_db->getScalar("SizeStep"); + if (rev_db->keyExists("CCInterval")) + cc_interval = rev_db->getScalar("CCInterval"); + if (rev_db->keyExists("gamma")) + gamma = rev_db->getScalar("gamma"); + if (rev_db->keyExists("StatSamples")) + stat_samples = rev_db->getScalar("StatSamples"); + + INSIST(size_step > 0.0, "REV.SizeStep must be positive"); + INSIST(cc_interval > 0, "REV.CCInterval must be positive"); + INSIST(stat_samples > 0, "REV.StatSamples must be positive"); + INSIST(std::isfinite(gamma) && gamma > 0.0, "REV.gamma must be positive and finite"); + auto simulation_params = mrt_db; + if (rev_db->keyExists("SimulationParams")) { + auto simulation_filename = + rev_db->getScalar("SimulationParams"); + INSIST(!simulation_filename.empty(), "REV.SimulationParams cannot be empty"); + const bool absolute = simulation_filename.front() == '/' || + simulation_filename.front() == '\\' || + (simulation_filename.size() > 1 && + simulation_filename[1] == ':'); + if (!absolute) { + auto separator = filename.find_last_of("/\\"); + if (separator != std::string::npos) + simulation_filename = filename.substr(0, separator + 1) + + simulation_filename; + } + + auto simulation_db = std::make_shared(simulation_filename); + INSIST(simulation_db->keyExists("MRT"), + "REV.SimulationParams must reference a database with an MRT section"); + simulation_params = simulation_db->getDatabase("MRT"); + } + if (simulation_params) { + tau = simulation_params->getWithDefault("tau", 1.0); + if (simulation_params->keyExists("F")) { + auto force = simulation_params->getVector("F"); + INSIST(force.size() == 3, "MRT.F must have three components"); + Fx = force[0]; + Fy = force[1]; + Fz = force[2]; + } + } + selection_flow_params = simulation_params; +} + +void ScaLBL_REV::SetDomain() { + standalone_ready = false; + stat_poro_rev = StatisticalREVResult{}; + stat_perm_rev = StatisticalREVResult{}; + if (!run_rev || (!geometrical && !hydraulic)) + return; + INSIST(db, "Call REV.ReadParams before REV.SetDomain"); + INSIST(comm.getSize() == 1, + "Standalone REV analysis currently requires one MPI process"); + + auto domain_db = db->getDatabase("Domain"); + Dm = std::make_shared(domain_db, comm); + Mask = Dm; + Nx = Dm->Nx; + Ny = Dm->Ny; + Nz = Dm->Nz; + N = Nx * Ny * Nz; + Lx = Dm->Lx; + Ly = Dm->Ly; + Lz = Dm->Lz; + rank = Dm->rank(); + nprocx = Dm->nprocx(); + nprocy = Dm->nprocy(); + nprocz = Dm->nprocz(); + + // Use the full domain when building the signed distance field. + std::fill(Dm->id.begin(), Dm->id.end(), 1); + Dm->CommInit(); + if (domain_db->keyExists("Filename")) { + Dm->Decomp(domain_db->getScalar("Filename")); + } else if (domain_db->keyExists("GridFile")) { + auto input_id = readMicroCT(*domain_db, comm); + INSIST(input_id.size(0) == size_t(Nx - 2) && + input_id.size(1) == size_t(Ny - 2) && + input_id.size(2) == size_t(Nz - 2), + "Input geometry does not match Domain.n"); + Array id_view; + id_view.viewRaw({size_t(Nx), size_t(Ny), size_t(Nz)}, Dm->id.data()); + fillHalo fill(comm, Dm->rank_info, + {Nx - 2, Ny - 2, Nz - 2}, {1, 1, 1}, 0, 1); + fill.copy(input_id, id_view); + } else { + Dm->ReadIDs(); + } + Np = Dm->PoreCount(); + + Distance.resize(Nx, Ny, Nz); + Array pore(Nx, Ny, Nz); + for (int i = 0; i < N; ++i) + pore(i) = Dm->id[i] > 0 ? 1 : 0; + CalcDist(Distance, pore, *Dm); + + Velocity_x.resize(Nx, Ny, Nz); + Velocity_y.resize(Nx, Ny, Nz); + Velocity_z.resize(Nx, Ny, Nz); + Velocity_x.fill(0.0); + Velocity_y.fill(0.0); + Velocity_z.fill(0.0); + if (hydraulic) + LoadVelocityFields(); + + average_pore_size = 0.0; + rev_x_poro = rev_x_perm = rev_x_tort = rev_x_surf = -1.0; + standalone_ready = true; +} + +void ScaLBL_REV::LoadVelocityFields() { + auto rev_db = db->getDatabase("REV"); + auto path = rev_db->getWithDefault("VelocityPath", "."); + std::string timestep; + if (rev_db->keyExists("VelocityTimestep")) { + timestep = rev_db->getScalar("VelocityTimestep"); + } else { + INSIST(exists(path + "/LBM.visit") || exists(path + "/summary.LBM"), + "Standalone hydraulic REV requires saved velocity output: " + "set REV.VelocityPath and REV.VelocityTimestep, " + "or set REV.Hydraulic = false"); + auto timesteps = IO::readTimesteps(path); + INSIST(!timesteps.empty(), "No saved velocity timesteps found"); + timestep = timesteps.back(); + } + INSIST(exists(path + "/" + timestep + "/LBM.summary"), + "Cannot find LBM.summary for REV velocity input"); + auto meshes = IO::getMeshList(path, timestep); + auto entry = std::find_if(meshes.begin(), meshes.end(), + [](const IO::MeshDatabase& mesh) { + return mesh.name == "domain"; + }); + INSIST(entry != meshes.end(), "Saved velocity output has no domain mesh"); + INSIST(entry->domains.size() == 1, + "Standalone REV requires saved output from one MPI domain"); + auto mesh = IO::getMesh(path, timestep, *entry, 0); + auto domain_mesh = std::dynamic_pointer_cast(mesh); + INSIST(domain_mesh && domain_mesh->nx == Nx - 2 && + domain_mesh->ny == Ny - 2 && domain_mesh->nz == Nz - 2, + "Saved velocity dimensions do not match the REV geometry"); + + fillHalo fill(comm, Dm->rank_info, + {Nx - 2, Ny - 2, Nz - 2}, {1, 1, 1}, 0, 1); + const char* names[] = {"Velocity_x", "Velocity_y", "Velocity_z"}; + DoubleArray* fields[] = {&Velocity_x, &Velocity_y, &Velocity_z}; + for (int i = 0; i < 3; ++i) { + auto variable = IO::getVariable(path, timestep, *entry, 0, names[i]); + INSIST(variable && variable->dim == 1 && + variable->type == IO::VariableType::VolumeVariable, + std::string("Missing scalar velocity field: ") + names[i]); + IO::reformatVariable(*mesh, *variable); + INSIST(variable->data.size(0) == size_t(Nx - 2) && + variable->data.size(1) == size_t(Ny - 2) && + variable->data.size(2) == size_t(Nz - 2), + "Saved velocity array dimensions do not match the REV geometry"); + fill.copy(variable->data, *fields[i]); + } +} + +void ScaLBL_REV::PoreSize() { + INSIST(standalone_ready, "Call REV.SetDomain before REV.PoreSize"); + PoreSize(*this); +} + +void ScaLBL_REV::DetRevAnalysis() { + if (!run_rev || (!geometrical && !hydraulic)) + return; + INSIST(standalone_ready, "Call REV.SetDomain before REV.DetRevAnalysis"); + DetRevAnalysis(*this); +} + +void ScaLBL_REV::StatRevAnalysis() { + if (!run_rev || (!geometrical && !hydraulic)) + return; + INSIST(standalone_ready, "Call REV.SetDomain before REV.StatRevAnalysis"); + if (average_pore_size == 0.0) + PoreSize(); + StatRevAnalysis(*this); +} + +void ScaLBL_REV::sub_domain_select() { + if (!run_rev || !select_subdomains || (!geometrical && !hydraulic)) + return; + INSIST(standalone_ready, "Call REV.SetDomain before REV.sub_domain_select"); + sub_domain_select(*this); +} + +struct det_convergence_metrics { + double RE; + double CC; +}; + +struct stat_convergence_metrics { + double CV, mean, max, min; +}; + +struct VoxelDist { + int x, y, z; + double dist; +}; + +det_convergence_metrics det_compute_convergence(const std::vector& data, + const std::vector& volumes, + int iter_step, int cc_interval) { + det_convergence_metrics m = {-1, -1}; + + if (iter_step < 1) + return m; + + double mean = 0.0; + double std_dev = 0.0; + int valid_count = 0; + int start_idx = (iter_step >= cc_interval - 1) ? (iter_step + 1 - cc_interval) : 0; + for (int l = start_idx; l <= iter_step; l++) { + if (data[l] != -1) { + mean += data[l]; + valid_count++; + } + } + + if (valid_count == 0) return m; + + mean /= valid_count; + + for (int l = start_idx; l <= iter_step; l++) { + if (data[l] != -1) { + std_dev += (data[l] - mean) * (data[l] - mean); + } + } + std_dev = sqrt(std_dev / valid_count); + + if (mean != 0) { + m.CC = std_dev / mean; + } + + if (data[iter_step] != -1 && data[iter_step - 1] != -1) { + double sum_prev_curr = data[iter_step] + data[iter_step - 1]; + if (sum_prev_curr != 0 && volumes[iter_step] > 0) { + m.RE = fabs((data[iter_step] - data[iter_step - 1]) / sum_prev_curr * + (volumes[iter_step - 1] / volumes[iter_step])); + } + } + + return m; +} + +stat_convergence_metrics stat_compute_convergence(const std::vector& data, int samples) { + stat_convergence_metrics res = {-1.0, 0, -1, 1e100}; + int valid_count = 0; + + for (int l = 0; l < samples; l++) { + if (data[l] != -1) { + res.mean += data[l]; + if (data[l] > res.max) res.max = data[l]; + if (data[l] < res.min) res.min = data[l]; + valid_count++; + } + } + + if (valid_count == 0) return res; + + res.mean /= valid_count; + + double std_dev = 0.0; + for (int l = 0; l < samples; l++) { + if (data[l] != -1) { + std_dev += (data[l] - res.mean) * (data[l] - res.mean); + } + } + std_dev = std::sqrt(std_dev / valid_count); + + if (res.mean != 0) res.CV = std_dev / std::abs(res.mean); + else if (std_dev == 0.0) res.CV = 0.0; + + return res; +} + +void ScaLBL_REV::PoreSizeImpl(const ModelData &model) { + int Nx = model.Nx, Ny = model.Ny, Nz = model.Nz; + double h = model.voxel_length; + + DoubleArray distance_updated; + distance_updated.resize(Nx, Ny, Nz); + + std::vector pore_voxels; + + for (int k = 0; k < Nz; k++) { + for (int j = 0; j < Ny; j++) { + for (int i = 0; i < Nx; i++) { + double val = model.Distance(i, j, k); + distance_updated(i, j, k) = val; + + if (i > 0 && i < Nx - 1 && j > 0 && j < Ny - 1 && k > 0 && k < Nz - 1) { + if (val > 0) { + pore_voxels.push_back({i, j, k, val}); + } + } + } + } + } + + std::sort(pore_voxels.begin(), pore_voxels.end(), [](const VoxelDist& a, const VoxelDist& b) { + return a.dist > b.dist; + }); + + for (const auto& voxel : pore_voxels) { + int i = voxel.x; + int j = voxel.y; + int k = voxel.z; + double current_max = voxel.dist; + + int start_x = std::max(1, static_cast(i - current_max - 1)); + int end_x = std::min(Nx - 1, static_cast(i + current_max + 2)); + int start_y = std::max(1, static_cast(j - current_max - 1)); + int end_y = std::min(Ny - 1, static_cast(j + current_max + 2)); + int start_z = std::max(1, static_cast(k - current_max - 1)); + int end_z = std::min(Nz - 1, static_cast(k + current_max + 2)); + + double r_sq = current_max * current_max; + + for (int n = start_z; n < end_z; n++) { + for (int m = start_y; m < end_y; m++) { + for (int l = start_x; l < end_x; l++) { + double dist_sq = (i - l) * (i - l) + (j - m) * (j - m) + (k - n) * (k - n); + + if (dist_sq <= r_sq) { + double value2 = distance_updated(l, m, n); + if (value2 > 0 && value2 < current_max) { + distance_updated(l, m, n) = current_max; + } + } + } + } + } + } + + double sum = 0.0; + int total_pore_voxels = 0; + + std::map psd_histogram; + + for (int k = 1; k < Nz - 1; k++) { + for (int j = 1; j < Ny - 1; j++) { + for (int i = 1; i < Nx - 1; i++) { + double val = distance_updated(i, j, k); + if (val > 0) { + sum += val; + total_pore_voxels++; + psd_histogram[val]++; + } + } + } + } + + average_pore_size = (total_pore_voxels > 0) ? (2.0 * sum / total_pore_voxels) : 0.0; + + std::ofstream outfile("pore_size_distribution.csv"); + if (outfile.is_open()) { + outfile << "radius_voxels diameter_um voxel_count volume_fraction\n"; + + for (const auto& pair : psd_histogram) { + double radius_vx = pair.first; + int count = pair.second; + + double diameter_um = 2.0 * radius_vx * h; + double volume_fraction = static_cast(count) / total_pore_voxels; + + outfile << radius_vx << " " + << diameter_um << " " + << count << " " + << volume_fraction << "\n"; + } + + outfile << "\nAverage_Pore_Size_um " << average_pore_size * h << "\n"; + outfile.close(); + } + + std::cout << "\n\n\nAverage Pore Size: " << average_pore_size * h << " micrometers\n\n\n" << std::endl; +} + +void ScaLBL_REV::DetRevAnalysisImpl(const ModelData &model) { + if (!run_rev) + return; + + PoreSizeImpl(model); + + double current_size_x_double = 2.0 * average_pore_size; + int x_size = static_cast(std::ceil(current_size_x_double)); //could be a problem if the geometry is in the yz plane. + + + int Nx = model.Nx, Ny = model.Ny, Nz = model.Nz; + int iter_step = 0; + double h = model.voxel_length; + double mu = (model.tau - 0.5) / 3.f; + + rev_x_poro = -1; + rev_x_perm = -1; + rev_x_tort = -1; + rev_x_surf = -1; + + int vector_size_x = 1000; + std::vector sub_volume(vector_size_x); + std::vector perm(vector_size_x); + std::vector poro(vector_size_x); + std::vector tort(vector_size_x); + std::vector surf(vector_size_x); + + std::ofstream geometric_log; + if (geometrical) { + geometric_log.open("det_rev_geometric.csv", std::ios::app); + if (geometric_log.tellp() == 0) { + geometric_log << "iter_step x_size_um x_size_vx " + << "porosity rev_poro_um rev_poro_vx RE_poro CC_poro " + << "surface rev_surf_um rev_surf_vx RE_surf CC_surf\n"; + } + } + + std::ofstream hydraulic_log; + if (hydraulic) { + hydraulic_log.open("det_rev_hydraulic.csv", std::ios::app); + if (hydraulic_log.tellp() == 0) { + hydraulic_log << "iter_step x_size_um x_size_vx " + << "permeability rev_perm_um rev_perm_vx RE_perm CC_perm " + << "tortuosity rev_tort_um rev_tort_vx RE_tort CC_tort\n"; + } + } + + + const int max_x_size = Nx - 2; + while (true) { + x_size = std::min(x_size, max_x_size); + + //(still have to test for 2d geometries) + int y_size = (Ny - 2) * x_size / (Nx - 2); + int z_size = (Nz - 2) * x_size / (Nx - 2); + + int start_x = (Nx - 2 - x_size) / 2 + 1; + int start_y = (Ny - 2 - y_size) / 2 + 1; + int start_z = (Nz - 2 - z_size) / 2 + 1; + + int end_x = (Nx - 2 + x_size) / 2 + 1; + int end_y = (Ny - 2 + y_size) / 2 + 1; + int end_z = (Nz - 2 + z_size) / 2 + 1; + + start_x = std::max(1, start_x); + end_x = std::min(Nx-2, end_x); + start_y = std::max(1, start_y); + end_y = std::min(Ny-2, end_y); + start_z = std::max(1, start_z); + end_z = std::min(Nz-2, end_z); + + // The common voxel volume cancels in V_(i-1) / V_i. + sub_volume[iter_step] = static_cast(end_x - start_x + 1) * + (end_y - start_y + 1) * (end_z - start_z + 1); + + double current_surf = 0; + double count = 0.0; + double vxb = 0.0, vyb = 0.0, vzb = 0.0; + double vax = 0.0, vay = 0.0, vaz = 0.0; + double v_tort=0.0; + + for (int k = start_z; k <= end_z; k++) { + for (int j = start_y; j <= end_y; j++) { + for (int i = start_x; i <= end_x; i++) { + if (model.Distance(i, j, k) > 0) { + if (hydraulic) { + vxb = model.Velocity_x(i, j, k); + vyb = model.Velocity_y(i, j, k); + vzb = model.Velocity_z(i, j, k); + vax += vxb; + vay += vyb; + vaz += vzb; + v_tort += sqrt(vxb*vxb + vyb*vyb + vzb*vzb); + } + count += 1.0; + } + if (geometrical) { + if ((model.Distance(i, j, k) * model.Distance(i, j, k + 1) < 0) && (k != end_z)) + current_surf++; + + if ((model.Distance(i, j, k) * model.Distance(i, j + 1, k) < 0) && (j != end_y)) + current_surf++; + + if ((model.Distance(i, j, k) * model.Distance(i + 1, j, k) < 0) && (i != end_x)) + current_surf++; + } + } + } + } + + if (hydraulic) { + if (count == 0) + { + vax = 0.0; + vay = 0.0; + vaz = 0.0; + } else { + vax /= count; + vay /= count; + vaz /= count;} + + double force_mag = sqrt(model.Fx * model.Fx + model.Fy * model.Fy + model.Fz * model.Fz); + double dir_x = model.Fx / force_mag; + double dir_y = model.Fy / force_mag; + double dir_z = model.Fz / force_mag; + if (force_mag == 0.0) { + dir_x = 0.0; + dir_y = 0.0; + dir_z = 1.0; + force_mag = 1.0; + } + + double flow_rate = (vax * dir_x + vay * dir_y + vaz * dir_z); + double current_poro = count / ((end_x - start_x + 1) * (end_y - start_y + 1) * (end_z - start_z + 1)); + double current_perm = 1013 * h * h * mu * current_poro * flow_rate / force_mag; + double current_tort = (vaz != 0) ? (v_tort / (vaz * count) - 1.0) : -1; + + perm[iter_step] = current_perm; + tort[iter_step] = current_tort; + + //Permeability deterministic approach + det_convergence_metrics perm_m = det_compute_convergence(perm, sub_volume, iter_step, cc_interval); + if (perm_m.RE >= 0 && perm_m.RE < gamma && perm_m.CC >= 0 && perm_m.CC < gamma && rev_x_perm == -1) { + rev_x_perm = x_size; + } + + //Tortuosity deterministic approach + det_convergence_metrics tort_m = det_compute_convergence(tort, sub_volume, iter_step, cc_interval); + if (tort_m.RE >= 0 && tort_m.RE < gamma && tort_m.CC >= 0 && tort_m.CC < gamma && rev_x_tort == -1) { + rev_x_tort = x_size; + } + + hydraulic_log << iter_step << " " << x_size * h << " " << x_size << " " + << perm[iter_step] << " " << rev_x_perm*h << " " << rev_x_perm << " " << perm_m.RE << " " << perm_m.CC << " " + << tort[iter_step]+1 << " " << rev_x_tort*h << " " << rev_x_tort << " " << tort_m.RE << " " << tort_m.CC << "\n"; + } + + if (geometrical) { + double vol = (end_x - start_x + 1) * (end_y - start_y + 1) * (end_z - start_z + 1); + current_surf = (current_surf / vol) / h; + double current_poro = count / vol; + + poro[iter_step] = current_poro; + surf[iter_step] = current_surf; + + //Porosity deterministic approach + det_convergence_metrics poro_m = det_compute_convergence(poro, sub_volume, iter_step, cc_interval); + if (poro_m.RE >= 0 && poro_m.RE < gamma && poro_m.CC >= 0 && poro_m.CC < gamma && rev_x_poro == -1) { + rev_x_poro = x_size; + } + + //Surface Deterministic approach + det_convergence_metrics surf_m = det_compute_convergence(surf, sub_volume, iter_step, cc_interval); + if (surf_m.RE >= 0 && surf_m.RE < gamma && surf_m.CC >= 0 && surf_m.CC < gamma && rev_x_surf == -1) { + rev_x_surf = x_size; + } + + geometric_log << iter_step << " " << x_size * h << " " << x_size << " " + << poro[iter_step] << " " << rev_x_poro*h << " " << rev_x_poro << " " << poro_m.RE << " " << poro_m.CC << " " + << surf[iter_step] << " " << rev_x_surf*h << " " << rev_x_surf << " " << surf_m.RE << " " << surf_m.CC << "\n"; + } + + // Stop after processing the full domain once. + if (x_size == max_x_size) { + break; + } + + double min_step = std::max(1.0, average_pore_size); + double proportional_step = current_size_x_double * size_step; + current_size_x_double += std::max(min_step, proportional_step); + x_size = static_cast(std::ceil(current_size_x_double)); + + iter_step++; + } + + if (geometrical) { + geometric_log.close(); + } + + if (hydraulic) { + hydraulic_log.close(); + } + + printf("\n\nEnded Deterministic Analysis \n\n"); +} + +void ScaLBL_REV::StatRevAnalysisImpl(const ModelData &model) { + if (!run_rev) + return; + + stat_poro_rev = StatisticalREVResult{}; + stat_perm_rev = StatisticalREVResult{}; + + std::mt19937 rng(0); + if (!seed_anchor) { + std::random_device rd; + rng.seed(rd()); + } + + int Nx = model.Nx, Ny = model.Ny, Nz = model.Nz; + + double h = model.voxel_length; + double mu = (model.tau - 0.5) / 3.0f; + + std::vector stat_poro(stat_samples); + std::vector stat_perm(stat_samples); + std::vector stat_tort(stat_samples); + std::vector stat_surf(stat_samples); + + std::ofstream geometric_log; + if (geometrical) { + geometric_log.open("stat_rev_geometric.csv", std::ios::app); + if (geometric_log.tellp() == 0) { + geometric_log << "x_size_um x_size_vx " + << "stat_poro stat_poro_mean stat_poro_max stat_poro_min " + << "stat_surf stat_surf_mean stat_surf_max stat_surf_min\n"; + } + } + + std::ofstream hydraulic_log; + if (hydraulic) { + hydraulic_log.open("stat_rev_hydraulic.csv", std::ios::app); + if (hydraulic_log.tellp() == 0) { + hydraulic_log << "x_size_um x_size_vx " + << "stat_perm stat_perm_mean stat_perm_max stat_perm_min " + << "stat_tort stat_tort_mean stat_tort_max stat_tort_min\n"; + } + } + + double current_size_x_double = 1e100; + if (rev_x_poro > 0) current_size_x_double = std::min(current_size_x_double, (double)rev_x_poro); + if (rev_x_perm > 0) current_size_x_double = std::min(current_size_x_double, (double)rev_x_perm); + if (rev_x_tort > 0) current_size_x_double = std::min(current_size_x_double, (double)rev_x_tort); + if (rev_x_surf > 0) current_size_x_double = std::min(current_size_x_double, (double)rev_x_surf); + if (current_size_x_double == 1e100) { + current_size_x_double = 2.0 * average_pore_size; + } + if (current_size_x_double == -1) current_size_x_double = 2.0 * average_pore_size; + int x_size = static_cast(std::ceil(current_size_x_double)); + + while (x_size > 0 && x_size < Nx - 2) { + + // Match the reported edge length: a cubic image uses L x L x L. + int sub_y = std::max(1, int(int64_t(x_size) * (Ny - 2) / (Nx - 2))); + int sub_z = std::max(1, int(int64_t(x_size) * (Nz - 2) / (Nx - 2))); + + int range_x = (Nx - 2) - x_size + 1; + int range_y = (Ny - 2) - sub_y + 1; + int range_z = (Nz - 2) - sub_z + 1; + + if (range_x < 1) range_x = 1; + if (range_y < 1) range_y = 1; + if (range_z < 1) range_z = 1; + + std::uniform_int_distribution dist_x(1, range_x); + std::uniform_int_distribution dist_y(1, range_y); + std::uniform_int_distribution dist_z(1, range_z); + + for (int l = 0; l < stat_samples; l++) { + int start_x = dist_x(rng); + int start_y = dist_y(rng); + int start_z = dist_z(rng); + + int end_x = std::min(Nx - 2, start_x + x_size - 1); + int end_y = std::min(Ny - 2, start_y + sub_y - 1); + int end_z = std::min(Nz - 2, start_z + sub_z - 1); + + double vax = 0.0, vay = 0.0, vaz = 0.0, count = 0.0, v_tort = 0.0, current_surf = 0.0; + + for (int k = start_z; k <= end_z; k++) { + for (int j = start_y; j <= end_y; j++) { + for (int i = start_x; i <= end_x; i++) { + double dist_val = model.Distance(i, j, k); + if (dist_val > 0) { + if (hydraulic) { + vax += model.Velocity_x(i, j, k); + vay += model.Velocity_y(i, j, k); + vaz += model.Velocity_z(i, j, k); + v_tort += std::sqrt(std::pow(model.Velocity_x(i, j, k), 2) + + std::pow(model.Velocity_y(i, j, k), 2) + + std::pow(model.Velocity_z(i, j, k), 2)); + } + count += 1.0; + } + + if (geometrical) { + if (k != end_z && (dist_val * model.Distance(i, j, k + 1) < 0)) current_surf++; + if (j != end_y && (dist_val * model.Distance(i, j + 1, k) < 0)) current_surf++; + if (i != end_x && (dist_val * model.Distance(i + 1, j, k) < 0)) current_surf++; + } + } + } + } + + double vol = (end_x - start_x + 1) * (end_y - start_y + 1) * (end_z - start_z + 1); + double current_poro = count / vol; + + if (hydraulic) { + if (count > 0) { + vax /= count; vay /= count; vaz /= count; + } + + double force_mag = std::sqrt(model.Fx * model.Fx + model.Fy * model.Fy + model.Fz * model.Fz); + double dir_x = (force_mag > 0) ? model.Fx / force_mag : 0.0; + double dir_y = (force_mag > 0) ? model.Fy / force_mag : 0.0; + double dir_z = (force_mag > 0) ? model.Fz / force_mag : 1.0; + if (force_mag == 0.0) force_mag = 1.0; + + double flow_rate = (vax * dir_x + vay * dir_y + vaz * dir_z); + double current_perm = 1013.0 * h * h * mu * current_poro * flow_rate / force_mag; + double current_tort = (vaz != 0) ? (v_tort / (vaz * count) - 1.0) : -1; + + stat_perm[l] = current_perm; + stat_tort[l] = current_tort; + } + + if (geometrical) { + current_surf /= vol; + current_surf /= h; + stat_poro[l] = current_poro; + stat_surf[l] = current_surf; + } + } + + if (hydraulic) { + stat_convergence_metrics k_stats = stat_compute_convergence(stat_perm, stat_samples); + stat_convergence_metrics t_stats = stat_compute_convergence(stat_tort, stat_samples); + const bool has_tortuosity_samples = + std::any_of(stat_tort.begin(), stat_tort.end(), [](double value) { + return value != -1.0; + }); + if (stat_perm_rev.size == 0 && k_stats.CV >= 0.0 && k_stats.CV < gamma && + std::isfinite(k_stats.mean)) { + stat_perm_rev = {x_size, {{x_size, sub_y, sub_z}}, k_stats.mean, k_stats.CV}; + // The existing statistical samples store T - 1, using -1 for + // undefined flow. Do not turn an empty sample set into T = 1. + if (std::isfinite(t_stats.mean) && + std::any_of(stat_tort.begin(), stat_tort.end(), [](double value) { + return value != -1.0 && std::isfinite(value); + })) + stat_perm_rev.tortuosity_mean = t_stats.mean + 1.0; + } + hydraulic_log << x_size * h << " " << x_size << " " + << k_stats.CV << " " << k_stats.mean << " " << k_stats.max << " " << k_stats.min << " " + << t_stats.CV << " " << (has_tortuosity_samples ? t_stats.mean + 1 : 0) << " " + << t_stats.max + 1 << " " << (has_tortuosity_samples ? t_stats.min + 1 : 0) << "\n"; + } + + if (geometrical) { + stat_convergence_metrics p_stats = stat_compute_convergence(stat_poro, stat_samples); + stat_convergence_metrics s_stats = stat_compute_convergence(stat_surf, stat_samples); + if (stat_poro_rev.size == 0 && p_stats.CV >= 0.0 && p_stats.CV < gamma && + std::isfinite(p_stats.mean)) { + stat_poro_rev = {x_size, {{x_size, sub_y, sub_z}}, p_stats.mean, p_stats.CV, + s_stats.mean}; + } + geometric_log << x_size * h << " " << x_size << " " + << p_stats.CV << " " << p_stats.mean << " " << p_stats.max << " " << p_stats.min << " " + << s_stats.CV << " " << s_stats.mean << " " << s_stats.max << " " << s_stats.min << "\n"; + } + + double min_step = std::max(1.0, average_pore_size); + double proportional_step = current_size_x_double * size_step; + current_size_x_double += std::max(min_step, proportional_step); + x_size = static_cast(std::ceil(current_size_x_double)); + } + + if (geometrical) + geometric_log.close(); + if (hydraulic) + hydraulic_log.close(); + printf("\n\nEnded Statistical Analysis\n\n"); +} diff --git a/analysis/rev_analysis.h b/analysis/rev_analysis.h new file mode 100644 index 00000000..2c24f7e7 --- /dev/null +++ b/analysis/rev_analysis.h @@ -0,0 +1,143 @@ + +#ifndef REV_ANALYSIS_H +#define REV_ANALYSIS_H +/* + Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University + Copyright Equnior ASA + + This file is part of the Open Porous Media project (OPM). + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ +/* + * Representative elementary volume analysis + */ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "common/ScaLBL.h" +#include "common/Communication.h" +#include "common/MPI.h" +#include "analysis/Minkowski.h" +#include "ProfilerApp.h" +#include "common/Domain.h" + + + +class ScaLBL_REV { +public: + ScaLBL_REV(int RANK, int NP, const Utilities::MPI &COMM); + ~ScaLBL_REV(); + + // Standalone analysis: read settings, load the domain, then run analysis. + // The standalone executable supplies false; explicit REV settings override it. + void ReadParams(string filename, bool hydraulic_default = true); + void SetDomain(); + void PoreSize(); + void DetRevAnalysis(); + void StatRevAnalysis(); + void sub_domain_select(); + + // Captured at the first statistical size with CV < REV.gamma. + // Keep these separate from the deterministic rev_x_* estimates. + struct StatisticalREVResult { + int size = 0; + std::array dimensions{{0, 0, 0}}; + double mean = 0.0; + double cv = -1.0; + // Specific surface measured in the same samples at the porosity REV. + double surface_mean = std::numeric_limits::quiet_NaN(); + // Displayed hydraulic tortuosity (not tortuosity - 1), at the permeability REV. + double tortuosity_mean = std::numeric_limits::quiet_NaN(); + }; + StatisticalREVResult stat_poro_rev, stat_perm_rev; + + // Models supply dimensions, domain, distance, velocity, tau and force. + template void PoreSize(const Model& model) { + PoreSizeImpl(GetModelData(model)); + } + template void DetRevAnalysis(const Model& model) { + if (run_rev && (geometrical || hydraulic)) + DetRevAnalysisImpl(GetModelData(model)); + } + template void StatRevAnalysis(const Model& model) { + if (run_rev && (geometrical || hydraulic)) + StatRevAnalysisImpl(GetModelData(model)); + } + template void sub_domain_select(const Model& model) { + if (run_rev && select_subdomains && (geometrical || hydraulic)) + SubDomainSelectImpl(GetModelData(model)); + } + + double average_pore_size; + double rev_x_poro, rev_x_perm, rev_x_tort, rev_x_surf; + double size_step; + int cc_interval, stat_samples; + double gamma; + bool geometrical = true; + bool hydraulic = true; + bool seed_anchor = true; + + DoubleArray Distance, Velocity_x, Velocity_y, Velocity_z; + double tau, Fx, Fy, Fz; + + int Nx, Ny, Nz, N, Np; + int rank, nprocx, nprocy, nprocz, nprocs; + double Lx, Ly, Lz; + + std::shared_ptr Dm; // this domain is for analysis + std::shared_ptr Mask; // this domain is for lbm + std::shared_ptr ScaLBL_Comm; + // input database + std::shared_ptr db; + +private: + // Borrow model arrays so generic calls do not copy simulation fields. + struct ModelData { + int Nx, Ny, Nz; + double voxel_length; + const DoubleArray& Distance; + const DoubleArray& Velocity_x; + const DoubleArray& Velocity_y; + const DoubleArray& Velocity_z; + double tau, Fx, Fy, Fz; + }; + + template static ModelData GetModelData(const Model& model) { + INSIST(model.Dm, "REV analysis requires an initialized domain"); + return {model.Nx, model.Ny, model.Nz, model.Dm->voxel_length, + model.Distance, model.Velocity_x, model.Velocity_y, + model.Velocity_z, model.tau, model.Fx, model.Fy, model.Fz}; + } + + void PoreSizeImpl(const ModelData& model); + void DetRevAnalysisImpl(const ModelData& model); + void StatRevAnalysisImpl(const ModelData& model); + void LoadVelocityFields(); + void SubDomainSelectImpl(const ModelData& model); + std::shared_ptr selection_flow_params; + Utilities::MPI comm; + bool run_rev = false; + bool standalone_ready = false; + bool select_subdomains = false; + bool save_subdomain_raw = false; + // Zero means use each property's retained statistical REV dimensions. + std::array subdomain_size{{0, 0, 0}}; +}; + +#endif // REV_ANALYSIS_H diff --git a/analysis/rev_subdomain.cpp b/analysis/rev_subdomain.cpp new file mode 100644 index 00000000..b342c00e --- /dev/null +++ b/analysis/rev_subdomain.cpp @@ -0,0 +1,552 @@ +#include "analysis/rev_analysis.h" +#include "analysis/rev_subdomain_flow.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { +using Shape = std::array; +constexpr double porosity_weight = 0.6; +constexpr double surface_weight = 0.4; +constexpr double permeability_weight = 0.6; +constexpr double tortuosity_weight = 0.4; + +double relativeDeviation(double local, double target) { + return target != 0.0 ? std::abs(local - target) / std::abs(target) + : (local == 0.0 ? 0.0 : std::numeric_limits::infinity()); +} + +struct Selection { + Shape start{{0, 0, 0}}; + double value = 0.0; + double deviation = std::numeric_limits::infinity(); + double secondary_value = 0.0; + double secondary_deviation = 0.0; + double score = std::numeric_limits::infinity(); + uint64_t candidates = 0; + bool found = false; + + void considerHydraulics(const Shape &position, double local, double target, + double local_tortuosity, double target_tortuosity) { + if (!std::isfinite(local)) + throw std::runtime_error("Subdomain selection received a nonfinite property"); + if (!std::isfinite(local_tortuosity)) { + ++candidates; // Undefined tortuosity is evaluated but cannot be ranked. + return; + } + const double relative = relativeDeviation(local, target); + const double tortuosity_relative = relativeDeviation(local_tortuosity, target_tortuosity); + update(position, local, relative, local_tortuosity, tortuosity_relative, + permeability_weight * relative + tortuosity_weight * tortuosity_relative); + } + + void considerGeometry(const Shape &position, double local, double target, + double local_surface, double target_surface) { + if (!std::isfinite(local) || !std::isfinite(local_surface)) + throw std::runtime_error("Subdomain selection received a nonfinite geometric property"); + const double relative = relativeDeviation(local, target); + const double surface_relative = relativeDeviation(local_surface, target_surface); + update(position, local, relative, local_surface, surface_relative, + porosity_weight * relative + surface_weight * surface_relative); + } + + void update(const Shape &position, double local, double relative, + double secondary, double secondary_relative, double candidate_score) { + ++candidates; + // Strict comparison preserves the first (z, y, x) position on ties. + if (!found || candidate_score < score) { + found = true; + start = position; + value = local; + deviation = relative; + secondary_value = secondary; + secondary_deviation = secondary_relative; + score = candidate_score; + } + } +}; + +std::string csvString(const std::string &value) { + std::string result = "\""; + for (char c : value) { + result += c; + if (c == '"') result += '"'; + } + return result + '"'; +} + +std::string rawFilename(const char *property, const Shape &size) { + return std::string("subdomain_") + property + '_' + std::to_string(size[0]) + 'x' + + std::to_string(size[1]) + 'x' + std::to_string(size[2]) + ".raw"; +} + +void protectSource(const std::string &filename, + const std::array &raw_files) { + // Check identity, including symlink and hard-link aliases, before any + // result or staging file can be opened for writing. + struct stat original; + if (::stat(filename.c_str(), &original) != 0) + return; // Missing sources are diagnosed when an export is needed. + const auto check_output = [&](const std::string &reserved) { + struct stat output; + if (::stat(reserved.c_str(), &output) == 0 && original.st_dev == output.st_dev && + original.st_ino == output.st_ino) + throw std::runtime_error("Selection output aliases the original image: " + reserved); + }; + for (const char *reserved : {"subdomain_selection.csv", "geometrical-rev-domain.db", + "hydraulic-rev-domain.db", "geometrical-rev-domain.db.tmp", + "hydraulic-rev-domain.db.tmp"}) { + check_output(reserved); + } + for (const auto &raw_file : raw_files) { + if (raw_file.empty()) continue; + check_output(raw_file); + check_output(raw_file + ".tmp"); + } +} + +// Domain::Decomp maps labels (and narrows 16-bit images). Export from the +// original file instead so selection never changes the source representation. +struct RawSource { + std::string filename, read_type; + Shape shape, offset{{0, 0, 0}}; + int bytes = 1; + std::ifstream input; + + RawSource(const Database &domain, const Shape &interior) { + if (domain.keyExists("Filename")) { + filename = domain.getScalar("Filename"); + read_type = domain.getWithDefault("ReadType", "8bit"); + if (read_type != "8bit" && read_type != "16bit") + throw std::runtime_error("Subdomain export supports raw 8bit/16bit geometry"); + bytes = read_type == "16bit" ? 2 : 1; + auto n = domain.getVector("N"); + if (n.size() != 3) + throw std::runtime_error("Subdomain export requires three Domain.N dimensions"); + std::copy(n.begin(), n.end(), shape.begin()); + if (domain.keyExists("offset")) { + auto origin = domain.getVector("offset"); + if (origin.size() != 3) + throw std::runtime_error("Domain.offset must have three coordinates"); + std::copy(origin.begin(), origin.end(), offset.begin()); + } + } else if (domain.keyExists("GridFile")) { + throw std::runtime_error("Subdomain export requires raw geometry or ID.00000; GridFile is unsupported"); + } else { + if (domain.getWithDefault("BC", 0) != 0) + throw std::runtime_error("Subdomain selection of ID.00000 requires BC = 0 to avoid implicit inlet/outlet layers"); + filename = "ID.00000"; + read_type = "signed8bit"; + for (int d = 0; d < 3; ++d) { + shape[d] = interior[d] + 2; + offset[d] = 1; + } + } + for (int d = 0; d < 3; ++d) + if (shape[d] <= 0 || offset[d] < 0 || + int64_t(offset[d]) + interior[d] > shape[d]) + throw std::runtime_error("Analyzed domain exceeds the original raw image bounds"); + input.open(filename, std::ios::binary | std::ios::ate); + const auto expected = std::streamoff(shape[0]) * shape[1] * shape[2] * bytes; + if (!input || input.tellg() != expected) + throw std::runtime_error("Raw source size does not match its dimensions and data type: " + filename); + } + + Array poreMask(const Database &domain, const Shape &interior) { + // Decomp changes inlet/outlet voxels after relabeling. Recover the + // pre-modification pore mask using exactly its signed-byte convention, + // including the narrowing of 16-bit labels before label mapping. + std::array labels; + for (size_t i = 0; i < labels.size(); ++i) + labels[i] = static_cast(i); + if (domain.keyExists("Filename")) { + const auto read = domain.getVector("ReadValues"); + const auto write = domain.getVector("WriteValues"); + if (read.size() != write.size()) + throw std::runtime_error("ReadValues and WriteValues must have equal lengths"); + std::array mapped{}; + for (size_t i = 0; i < read.size(); ++i) { + const auto index = static_cast(read[i]); + if (!mapped[index]) { + labels[index] = static_cast(write[i]); + mapped[index] = true; + } + } + } + Array pores(interior[0], interior[1], interior[2]); + std::vector row(size_t(interior[0]) * bytes); + for (int z = 0; z < interior[2]; ++z) { + for (int y = 0; y < interior[1]; ++y) { + const auto index = (std::streamoff(z + offset[2]) * shape[1] + + y + offset[1]) * shape[0] + offset[0]; + input.seekg(index * bytes); + input.read(row.data(), std::streamsize(row.size())); + if (!input) + throw std::runtime_error("Failed to read original rock pore labels"); + for (int x = 0; x < interior[0]; ++x) { + unsigned char label; + if (bytes == 2) { + uint16_t original; + std::memcpy(&original, row.data() + size_t(x) * 2, 2); + label = static_cast(original); + } else { + label = static_cast(row[x]); + } + pores(x, y, z) = labels[label] > 0 ? 1 : 0; + } + } + } + return pores; + } + + void save(const std::string &output, const Shape &start, const Shape &size) { + const std::string temporary = output + ".tmp"; + std::ofstream out(temporary, std::ios::binary | std::ios::trunc); + if (!out) throw std::runtime_error("Cannot write selected raw image: " + temporary); + std::vector row(size_t(size[0]) * bytes); + for (int z = 0; z < size[2]; ++z) { + for (int y = 0; y < size[1]; ++y) { + auto index = (std::streamoff(start[2] + offset[2] + z) * shape[1] + + start[1] + offset[1] + y) * shape[0] + start[0] + offset[0]; + input.seekg(index * bytes); + input.read(row.data(), std::streamsize(row.size())); + out.write(row.data(), std::streamsize(row.size())); + if (!input || !out) + throw std::runtime_error("Failed to copy original subdomain bytes"); + } + } + out.close(); + if (!out || std::rename(temporary.c_str(), output.c_str()) != 0) + throw std::runtime_error("Cannot finish selected raw image: " + output); + } + + void saveDomain(const std::string &output, const Database &original_domain, + const Shape &start, const Shape &size, double voxel_length) const { + auto selected = original_domain.cloneDatabase(); + char *resolved = ::realpath(filename.c_str(), nullptr); + if (!resolved) + throw std::runtime_error("Cannot resolve original geometry path: " + filename); + const std::string absolute_filename(resolved); + std::free(resolved); + selected->putScalar("Filename", absolute_filename); + selected->putScalar("ReadType", bytes == 2 ? "16bit" : "8bit"); + selected->putVector("N", {shape[0], shape[1], shape[2]}); + selected->putVector("n", {size[0], size[1], size[2]}); + selected->putVector("offset", {offset[0] + start[0], offset[1] + start[1], + offset[2] + start[2]}); + selected->putVector("nproc", {1, 1, 1}); + selected->putScalar("voxel_length", voxel_length); + if (selected->keyExists("L")) + selected->putVector("L", {size[0] * voxel_length, size[1] * voxel_length, + size[2] * voxel_length}); + selected->putScalar("BC", 0); + selected->putVector("InletLayers", {0, 0, 0}); + selected->putVector("OutletLayers", {0, 0, 0}); + if (!original_domain.keyExists("Filename")) { + // ID.00000 already stores mapped signed-byte labels. Reading it as + // a raw full image must preserve that mapping and include its halo + // in N and offset; unmatched labels also pass through unchanged. + selected->putVector("ReadValues", {0, 1, 2}); + selected->putVector("WriteValues", {0, 1, 2}); + } + Database result; + result.putDatabase("Domain", selected); + const std::string temporary = output + ".tmp"; + std::ofstream out(temporary, std::ios::trunc); + if (!out) + throw std::runtime_error("Cannot write selected Domain database: " + temporary); + out << std::setprecision(std::numeric_limits::max_digits10); + result.print(out); + out.close(); + if (!out || std::rename(temporary.c_str(), output.c_str()) != 0) + throw std::runtime_error("Cannot finish selected Domain database: " + output); + } +}; +} // namespace + +void ScaLBL_REV::SubDomainSelectImpl(const ModelData &model) { + if (comm.getSize() != 1) + throw std::runtime_error("Subdomain selection requires one MPI process"); + const Shape image{{model.Nx - 2, model.Ny - 2, model.Nz - 2}}; + const double h = model.voxel_length; + if (!std::isfinite(h) || h <= 0.0) + throw std::runtime_error("Subdomain selection requires a positive finite voxel length"); + auto domain = db->getDatabase("Domain"); + const auto configured_size = domain->getVector("n"); + if (configured_size.size() != 3 || + !std::equal(configured_size.begin(), configured_size.end(), image.begin())) + throw std::runtime_error("Selection model dimensions do not match Domain.n"); + if (subdomain_size[0] > 0) + for (int d = 0; d < 3; ++d) + if (subdomain_size[d] > image[d]) + throw std::runtime_error("REV.subdomain_size exceeds the analyzed Domain.n dimensions"); + const auto selection_dimensions = [&](const StatisticalREVResult &stat) { + return subdomain_size[0] > 0 ? subdomain_size : stat.dimensions; + }; + bool synthetic_layers = false; + for (const char *key : {"InletLayers", "OutletLayers"}) { + if (domain->keyExists(key)) { + for (int layers : domain->getVector(key)) + synthetic_layers = synthetic_layers || layers != 0; + } + } + std::array raw_files; + if (save_subdomain_raw) { + if (geometrical && stat_poro_rev.size > 0) + raw_files[0] = rawFilename("porosity", selection_dimensions(stat_poro_rev)); + if (hydraulic && stat_perm_rev.size > 0) + raw_files[1] = rawFilename("permeability", selection_dimensions(stat_perm_rev)); + } + std::unique_ptr source; + protectSource(domain->getWithDefault("Filename", "ID.00000"), raw_files); + if ((geometrical && stat_poro_rev.size > 0) || + (hydraulic && stat_perm_rev.size > 0)) + source.reset(new RawSource(*domain, image)); + + Array original_pores; + if (synthetic_layers && source) { + original_pores = source->poreMask(*domain, image); + std::cout << "Subdomain selection uses original rock labels before inlet/outlet modifications; " + "statistical REV targets are retained." << std::endl; + } + const auto is_pore = [&](int x, int y, int z) { + return synthetic_layers ? original_pores(x, y, z) > 0 + : model.Distance(x + 1, y + 1, z + 1) > 0.0; + }; + + std::ofstream summary("subdomain_selection.csv", std::ios::trunc); + if (!summary) throw std::runtime_error("Cannot write subdomain_selection.csv"); + summary << std::setprecision(std::numeric_limits::max_digits10) + << "property,status,rev_size_vx,rev_size_um,nx,ny,nz,lx_um,ly_um,lz_um," + "x0,y0,z0,x1,y1,z1,statistical_mean,selected_value,relative_deviation," + "candidate_count,raw_file,source_file,read_type,bytes_per_voxel," + "source_nx,source_ny,source_nz,source_x0,source_y0,source_z0," + "source_x1,source_y1,source_z1,domain_file,statistical_mean_surface," + "selected_surface,relative_deviation_surface,selection_score,porosity_weight,surface_weight," + "statistical_mean_tortuosity,selected_tortuosity,relative_deviation_tortuosity," + "permeability_weight,tortuosity_weight\n"; + + const auto write_status = [&](const char *name, const char *status, + const StatisticalREVResult *stat, uint64_t candidates) { + summary << name << ',' << status; + const Shape size = stat ? selection_dimensions(*stat) : Shape{{0, 0, 0}}; + for (int column = 2; column < 45; ++column) { + summary << ','; + if (column == 19) summary << candidates; + if (!stat) continue; + if (column == 2) summary << stat->size; + else if (column == 3) summary << stat->size * h; + else if (column >= 4 && column <= 6) summary << size[column - 4]; + else if (column >= 7 && column <= 9) summary << size[column - 7] * h; + else if (column == 16) summary << stat->mean; + else if (column == 40 && std::isfinite(stat->tortuosity_mean)) + summary << stat->tortuosity_mean; + } + summary << '\n'; + summary.flush(); + if (!summary) throw std::runtime_error("Failed to write subdomain_selection.csv"); + }; + + bool device_ready = false; + for (int property = 0; property < 2; ++property) { + const bool porosity = property == 0; + const char *name = porosity ? "porosity" : "permeability"; + const bool enabled = porosity ? geometrical : hydraulic; + const auto &stat = porosity ? stat_poro_rev : stat_perm_rev; + if (!enabled || stat.size == 0) { + write_status(name, enabled ? "not_converged" : "disabled", nullptr, 0); + continue; + } + const Shape size = selection_dimensions(stat); + for (int d = 0; d < 3; ++d) + if (size[d] < 1 || size[d] > image[d]) + throw std::runtime_error("Statistical REV dimensions exceed the input image"); + if (!std::isfinite(stat.mean)) + throw std::runtime_error("Statistical REV mean must be finite"); + if (!porosity && !std::isfinite(stat.tortuosity_mean)) { + write_status(name, "unavailable_tortuosity", &stat, 0); + std::cout << "Hydraulic selection has no valid statistical tortuosity target at the permeability REV size." + << std::endl; + continue; + } + const double volume = double(size[0]) * size[1] * size[2]; + const uint64_t expected = uint64_t(image[0] - size[0] + 1) * + (image[1] - size[1] + 1) * (image[2] - size[2] + 1); + std::cout << "Selecting " << name << " subdomain: " << expected + << " windows of " << size[0] << 'x' << size[1] << 'x' << size[2] << std::endl; + Selection best; + if (porosity) { + if (!std::isfinite(stat.surface_mean) || stat.surface_mean < 0.0) + throw std::runtime_error("Geometric selection requires the statistical surface mean at the porosity REV size"); + // Four sliding z slabs and summed-area tables count pore voxels + // and internal x/y/z interface faces, with O(Nx*Ny) memory. + const size_t pitch = size_t(image[0]) + 1; + std::array, 4> slab; + std::array, 4> prefix; + for (int channel = 0; channel < 4; ++channel) { + slab[channel].resize(size_t(image[0]) * image[1], 0); + prefix[channel].resize(pitch * (size_t(image[1]) + 1), 0); + } + const auto interface_face = [&](int axis, int x, int y, int z) { + Shape next{{x, y, z}}; + if (++next[axis] >= image[axis]) return false; + if (synthetic_layers) + return is_pore(x, y, z) != is_pore(next[0], next[1], next[2]); + const double a = model.Distance(x + 1, y + 1, z + 1); + const double b = model.Distance(next[0] + 1, next[1] + 1, next[2] + 1); + // Match statistical analysis: zero distance is not an interface. + return (a < 0.0 && b > 0.0) || (a > 0.0 && b < 0.0); + }; + auto add_plane = [&](int z, int sign) { + for (int y = 0; y < image[1]; ++y) + for (int x = 0; x < image[0]; ++x) { + const size_t index = size_t(y) * image[0] + x; + slab[0][index] += sign * is_pore(x, y, z); + slab[1][index] += sign * interface_face(0, x, y, z); + slab[2][index] += sign * interface_face(1, x, y, z); + } + }; + auto add_z_faces = [&](int z, int sign) { + for (int y = 0; y < image[1]; ++y) + for (int x = 0; x < image[0]; ++x) + slab[3][size_t(y) * image[0] + x] += sign * interface_face(2, x, y, z); + }; + for (int z = 0; z < size[2]; ++z) add_plane(z, 1); + for (int z = 0; z < size[2] - 1; ++z) add_z_faces(z, 1); + const auto rectangle_sum = [&](int channel, int x, int y, int width, int height) { + const auto &table = prefix[channel]; + return (table[size_t(y + height) * pitch + x + width] - + table[size_t(y) * pitch + x + width]) - + (table[size_t(y + height) * pitch + x] - + table[size_t(y) * pitch + x]); + }; + for (int z = 0; z <= image[2] - size[2]; ++z) { + if (z > 0) { + add_plane(z - 1, -1); + add_plane(z + size[2] - 1, 1); + if (size[2] > 1) { + add_z_faces(z - 1, -1); + add_z_faces(z + size[2] - 2, 1); + } + } + for (int channel = 0; channel < 4; ++channel) { + for (int y = 0; y < image[1]; ++y) { + uint64_t row = 0; + for (int x = 0; x < image[0]; ++x) { + row += slab[channel][size_t(y) * image[0] + x]; + prefix[channel][size_t(y + 1) * pitch + x + 1] = + prefix[channel][size_t(y) * pitch + x + 1] + row; + } + } + } + for (int y = 0; y <= image[1] - size[1]; ++y) { + for (int x = 0; x <= image[0] - size[0]; ++x) { + const uint64_t count = rectangle_sum(0, x, y, size[0], size[1]); + const uint64_t faces = rectangle_sum(1, x, y, size[0] - 1, size[1]) + + rectangle_sum(2, x, y, size[0], size[1] - 1) + + rectangle_sum(3, x, y, size[0], size[1]); + best.considerGeometry({{x, y, z}}, count / volume, stat.mean, + (faces / volume) / h, stat.surface_mean); + } + } + } + } else { + int bc = domain->getWithDefault("BC", 0); + if (selection_flow_params) { + if (selection_flow_params->keyExists("BC")) + bc = selection_flow_params->getScalar("BC"); + else if (selection_flow_params->keyExists("BoundaryCondition")) + bc = selection_flow_params->getScalar("BoundaryCondition"); + } + if (bc != 0) + throw std::runtime_error("Permeability selection requires periodic body-force flow (BC = 0)"); + const int max_steps = selection_flow_params + ? selection_flow_params->getWithDefault("timestepMax", 100000) : 100000; + const int interval = selection_flow_params + ? selection_flow_params->getWithDefault("analysis_interval", 1000) : 1000; + const double tolerance = selection_flow_params + ? selection_flow_params->getWithDefault("tolerance", 1e-8) : 1e-8; + if (!device_ready) { + ScaLBL_SetDevice(rank); + device_ready = true; + } + Array window(size[0], size[1], size[2]); + for (int z = 0; z <= image[2] - size[2]; ++z) { + for (int y = 0; y <= image[1] - size[1]; ++y) { + for (int x = 0; x <= image[0] - size[0]; ++x) { + for (int k = 0; k < size[2]; ++k) + for (int j = 0; j < size[1]; ++j) + for (int i = 0; i < size[0]; ++i) + window(i, j, k) = is_pore(x + i, y + j, z + k) ? 1 : 0; + rev_selection::HydraulicProperties properties; + try { + properties = rev_selection::simulateHydraulicProperties( + window, h, model.tau, {{model.Fx, model.Fy, model.Fz}}, + max_steps, interval, tolerance, comm); + } catch (const std::exception &error) { + std::ostringstream message; + message << "Permeability candidate (" << x << ',' << y << ',' << z + << "): " << error.what(); + throw std::runtime_error(message.str()); + } + best.considerHydraulics({{x, y, z}}, properties.permeability, stat.mean, + properties.tortuosity, stat.tortuosity_mean); + } + } + std::cout << "Permeability selection: " << best.candidates << '/' << expected + << " windows evaluated" << std::endl; + } + } + if (best.candidates != expected) + throw std::runtime_error("Subdomain selection did not evaluate every valid window"); + if (!best.found) { + write_status(name, "no_valid_tortuosity", &stat, best.candidates); + std::cout << "Hydraulic selection found no candidate with defined tortuosity after " + << best.candidates << " windows." << std::endl; + continue; + } + const std::string &raw_file = raw_files[property]; + if (!raw_file.empty()) + source->save(raw_file, best.start, size); + const std::string domain_file = porosity ? "geometrical-rev-domain.db" : "hydraulic-rev-domain.db"; + source->saveDomain(domain_file, *domain, best.start, size, h); + const double error = best.deviation; + summary << name << ",selected," << stat.size << ',' << stat.size * h; + for (int d : size) summary << ',' << d; + for (int d : size) summary << ',' << d * h; + for (int d : best.start) summary << ',' << d; + for (int d = 0; d < 3; ++d) summary << ',' << best.start[d] + size[d]; + summary << ',' << stat.mean << ',' << best.value << ',' << error << ',' << best.candidates + << ',' << csvString(raw_file) << ',' << csvString(source->filename) + << ',' << csvString(source->read_type) << ',' << source->bytes; + for (int d : source->shape) summary << ',' << d; + for (int d = 0; d < 3; ++d) summary << ',' << best.start[d] + source->offset[d]; + for (int d = 0; d < 3; ++d) summary << ',' << best.start[d] + source->offset[d] + size[d]; + summary << ',' << csvString(domain_file); + if (porosity) + summary << ',' << stat.surface_mean << ',' << best.secondary_value << ',' << best.secondary_deviation + << ',' << best.score << ',' << porosity_weight << ',' << surface_weight << ",,,,,"; + else + summary << ",,,," << best.score << ",,," << stat.tortuosity_mean << ',' << best.secondary_value + << ',' << best.secondary_deviation << ',' << permeability_weight << ',' << tortuosity_weight; + summary << '\n'; + summary.flush(); + if (!summary) throw std::runtime_error("Failed to write subdomain_selection.csv"); + std::cout << "Selected " << name << " window at (" << best.start[0] << ',' + << best.start[1] << ',' << best.start[2] << "), value=" << best.value + << ", relative deviation=" << error; + if (porosity) + std::cout << ", specific surface=" << best.secondary_value; + else + std::cout << ", tortuosity=" << best.secondary_value; + std::cout << ", selection score=" << best.score; + std::cout << std::endl; + } +} diff --git a/analysis/rev_subdomain_flow.cpp b/analysis/rev_subdomain_flow.cpp new file mode 100644 index 00000000..93384c43 --- /dev/null +++ b/analysis/rev_subdomain_flow.cpp @@ -0,0 +1,249 @@ +#include "analysis/rev_subdomain_flow.h" + +#include "common/Database.h" +#include "common/Domain.h" +#include "common/ScaLBL.h" + +#include +#include +#include +#include +#include +#include + +namespace rev_selection { +namespace { + +// Simulation storage must be released after every window, including failures. +template class DeviceBuffer { +public: + explicit DeviceBuffer(size_t count) { + void *storage = nullptr; + ScaLBL_AllocateDeviceMemory(&storage, count * sizeof(T)); + data_ = static_cast(storage); + if (!data_) + throw std::runtime_error("Cannot allocate REV window flow storage"); + } + ~DeviceBuffer() { ScaLBL_FreeDeviceMemory(data_); } + DeviceBuffer(const DeviceBuffer &) = delete; + DeviceBuffer &operator=(const DeviceBuffer &) = delete; + T *data() const { return data_; } + +private: + T *data_ = nullptr; +}; + +} // namespace + +HydraulicProperties simulateHydraulicProperties( + const Array &labels, double voxel_length, double tau, + std::array force, int max_steps, int interval, double tolerance, + const Utilities::MPI &comm) { + if (comm.getSize() != 1) + throw std::runtime_error("REV window flow requires one MPI process"); + if (!std::isfinite(voxel_length) || voxel_length <= 0.0 || + !std::isfinite(tau) || tau <= 0.5 || + !std::isfinite(tolerance) || tolerance <= 0.0 || tolerance >= 1.0 || + max_steps < 2 || interval < 1) + throw std::runtime_error("Invalid REV window flow simulation settings"); + for (double component : force) + if (!std::isfinite(component)) + throw std::runtime_error("REV window flow force must be finite"); + const double force_magnitude = + std::hypot(std::hypot(force[0], force[1]), force[2]); + if (!std::isfinite(force_magnitude) || force_magnitude <= 0.0) + throw std::runtime_error("REV window flow requires a nonzero body force"); + + const size_t limit = size_t(std::numeric_limits::max()); + size_t padded_volume = 1; + std::array size; + for (int axis = 0; axis < 3; ++axis) { + const size_t length = labels.size(axis); + if (length == 0 || length > limit - 2 || + padded_volume > limit / (length + 2)) + throw std::runtime_error("REV window exceeds flow solver size limits"); + padded_volume *= length + 2; + size[axis] = static_cast(length); + } + if (labels.ndim() > 3 || labels.length() > limit / 19 - 32) + throw std::runtime_error("REV window exceeds D3Q19 indexing limits"); + + int pore_count = 0; + for (size_t n = 0; n < labels.length(); ++n) + pore_count += labels(n) > 0 ? 1 : 0; + if (pore_count == 0) + return {0.0, std::numeric_limits::quiet_NaN()}; + + auto domain_db = std::make_shared(); + domain_db->putVector("n", {size[0], size[1], size[2]}); + domain_db->putVector("nproc", {1, 1, 1}); + domain_db->putScalar("BC", 0); + domain_db->putScalar("voxel_length", voxel_length); + auto domain = std::make_shared(domain_db, comm); + + // Match periodic MRT boundaries, including the geometry halo. No labels in + // the source window are changed and no inlet/outlet reservoirs are added. + const int nx = domain->Nx, ny = domain->Ny, nz = domain->Nz; + for (int k = 0; k < nz; ++k) { + const int z = (k + size[2] - 1) % size[2]; + for (int j = 0; j < ny; ++j) { + const int y = (j + size[1] - 1) % size[1]; + for (int i = 0; i < nx; ++i) { + const int x = (i + size[0] - 1) % size[0]; + domain->id[(k * ny + j) * nx + i] = labels(x, y, z); + } + } + } + domain->CommInit(); + auto communicator = std::make_shared(domain); + IntArray map(nx, ny, nz); + const int allocation_count = (pore_count / 16 + 2) * 16; + std::vector neighbors(size_t(18) * allocation_count, 0); + const int np = communicator->MemoryOptimizedLayoutAA( + map, neighbors.data(), domain->id.data(), pore_count, 1); + if (np > allocation_count) + throw std::runtime_error("REV window flow layout exceeds allocation"); + + DeviceBuffer device_neighbors(size_t(18) * np); + DeviceBuffer distributions(size_t(19) * np); + DeviceBuffer velocity(size_t(3) * np); + ScaLBL_CopyToDevice(device_neighbors.data(), neighbors.data(), + size_t(18) * np * sizeof(int)); + ScaLBL_D3Q19_Init(distributions.data(), np); + const double relaxation_a = 1.0 / tau; + const double relaxation_b = + 8.0 * (2.0 - relaxation_a) / (8.0 - relaxation_a); + // The momentum estimator reads post-collision j/rho - F/2. Initialize + // j=F/2 so the physical initial velocity is zero, including sealed pores + // whose bounce-back otherwise preserves an artificial odd/even oscillation. + ScaLBL_D3Q19_AAeven_MRT( + distributions.data(), 0, np, np, relaxation_a, relaxation_b, + 0.5 * force[0], 0.5 * force[1], 0.5 * force[2]); + ScaLBL_DeviceBarrier(); + + std::vector host_velocity(size_t(3) * np); + const double viscosity = (tau - 0.5) / 3.0; + const double porosity = double(pore_count) / double(labels.length()); + bool have_previous = false; + double previous_flow = 0.0, previous_speed = 0.0, previous_z = 0.0; + double previous_tortuosity = std::numeric_limits::quiet_NaN(); + int stable_checks = 0, last_check = 0; + + // Use the same odd/even AA updates and second-order momentum estimator as + // ScaLBL_MRTModel::Run. Check only after complete pairs of timesteps. + const int final_step = max_steps - max_steps % 2; + for (int step = 0; step < final_step;) { + communicator->SendD3Q19AA(distributions.data()); + ScaLBL_D3Q19_AAodd_MRT( + device_neighbors.data(), distributions.data(), + communicator->FirstInterior(), communicator->LastInterior(), np, + relaxation_a, relaxation_b, force[0], force[1], force[2]); + communicator->RecvD3Q19AA(distributions.data()); + ScaLBL_D3Q19_AAodd_MRT( + device_neighbors.data(), distributions.data(), 0, + communicator->LastExterior(), np, relaxation_a, relaxation_b, + force[0], force[1], force[2]); + ScaLBL_DeviceBarrier(); + + communicator->SendD3Q19AA(distributions.data()); + ScaLBL_D3Q19_AAeven_MRT( + distributions.data(), communicator->FirstInterior(), + communicator->LastInterior(), np, relaxation_a, relaxation_b, + force[0], force[1], force[2]); + communicator->RecvD3Q19AA(distributions.data()); + ScaLBL_D3Q19_AAeven_MRT( + distributions.data(), 0, communicator->LastExterior(), np, + relaxation_a, relaxation_b, force[0], force[1], force[2]); + ScaLBL_DeviceBarrier(); + step += 2; + if (step - last_check < interval && step != final_step) + continue; + last_check = step; + + ScaLBL_D3Q19_Momentum_2nd_order( + distributions.data(), velocity.data(), np, + force[0], force[1], force[2]); + ScaLBL_DeviceBarrier(); + ScaLBL_CopyToHost(host_velocity.data(), velocity.data(), + size_t(3) * np * sizeof(double)); + double projected_sum = 0.0, speed_sum = 0.0, z_sum = 0.0; + for (int k = 1; k < nz - 1; ++k) + for (int j = 1; j < ny - 1; ++j) + for (int i = 1; i < nx - 1; ++i) { + const int n = map(i, j, k); + if (n < 0) + continue; + std::array local_velocity; + for (int axis = 0; axis < 3; ++axis) { + const double component = host_velocity[size_t(axis) * np + n]; + if (!std::isfinite(component)) + throw std::runtime_error("Nonfinite REV window flow velocity"); + projected_sum += component * force[axis] / force_magnitude; + local_velocity[axis] = component; + } + speed_sum += std::hypot(std::hypot(local_velocity[0], local_velocity[1]), + local_velocity[2]); + z_sum += local_velocity[2]; + } + const double flow = projected_sum / pore_count; + const double mean_speed = speed_sum / pore_count; + const double mean_z = z_sum / pore_count; + if (!std::isfinite(flow) || !std::isfinite(mean_speed) || !std::isfinite(mean_z)) + throw std::runtime_error("Nonfinite REV window flow measurement"); + // Recover the same dimensionless tortuosity as statistical REV from + // this velocity field. Exclude unresolved z-flow residuals rather than + // interpreting roundoff in sealed or transverse-flow windows as flow. + const double zero_flow = tolerance * std::max(force_magnitude, mean_speed); + double tortuosity = std::numeric_limits::quiet_NaN(); + if (std::abs(mean_z) > zero_flow && mean_speed > zero_flow) { + tortuosity = mean_speed / mean_z; + if (!std::isfinite(tortuosity)) + throw std::runtime_error("Nonfinite REV window tortuosity"); + } + + // A force-sized floor permits sealed, zero-flow pore space to converge + // without division by zero. Require two consecutive stable intervals. + const auto stable = [&](double value, double previous) { + const double scale = std::max(force_magnitude, + std::max(std::abs(value), std::abs(previous))); + return std::abs(value - previous) <= tolerance * scale; + }; + const bool stable_tortuosity = + (std::isnan(tortuosity) && std::isnan(previous_tortuosity)) || + (std::isfinite(tortuosity) && std::isfinite(previous_tortuosity) && + std::abs(tortuosity - previous_tortuosity) <= tolerance * + std::max(1.0, std::max(std::abs(tortuosity), std::abs(previous_tortuosity)))); + if (have_previous && stable(flow, previous_flow) && + stable(mean_speed, previous_speed) && stable(mean_z, previous_z) && stable_tortuosity) + ++stable_checks; + else + stable_checks = 0; + if (stable_checks >= 2) { + if (flow < -tolerance * force_magnitude) + throw std::runtime_error("REV window flow opposes its driving force"); + const double permeability = 1013.0 * voxel_length * voxel_length * + viscosity * porosity * + std::max(0.0, flow) / force_magnitude; + if (!std::isfinite(permeability)) + throw std::runtime_error("Nonfinite REV window permeability"); + return {permeability, tortuosity}; + } + previous_flow = flow; + previous_speed = mean_speed; + previous_z = mean_z; + previous_tortuosity = tortuosity; + have_previous = true; + } + throw std::runtime_error("REV window permeability simulation did not converge within timestepMax"); +} + +double simulatePermeability(const Array &labels, + double voxel_length, double tau, + std::array force, int max_steps, + int interval, double tolerance, + const Utilities::MPI &comm) { + return simulateHydraulicProperties(labels, voxel_length, tau, force, + max_steps, interval, tolerance, comm).permeability; +} + +} // namespace rev_selection diff --git a/analysis/rev_subdomain_flow.h b/analysis/rev_subdomain_flow.h new file mode 100644 index 00000000..8882a4f3 --- /dev/null +++ b/analysis/rev_subdomain_flow.h @@ -0,0 +1,37 @@ +#ifndef LBPM_REV_SUBDOMAIN_FLOW_H +#define LBPM_REV_SUBDOMAIN_FLOW_H + +#include "common/Array.h" +#include "common/MPI.h" + +#include + +namespace rev_selection { + +struct HydraulicProperties { + double permeability; + double tortuosity; +}; + +// Run a fresh, periodic, body-force-driven D3Q19 MRT solve on a window. +// Labels contain the interior only, in x-fast order; positive labels are pores. +// The returned permeability uses the same mDarcy conversion as ScaLBL_MRTModel. +// Tortuosity matches the displayed statistical REV quantity sum(|u|)/sum(u_z). +// It is NaN when the window has no resolved net z flow, including solid and +// sealed pore space. Its sign follows the existing z-velocity convention. +// A nonfinite or unconverged simulation raises std::runtime_error. +HydraulicProperties simulateHydraulicProperties( + const Array &labels, double voxel_length, double tau, + std::array force, int max_steps, int interval, double tolerance, + const Utilities::MPI &comm); + +// Compatibility interface for callers that only need permeability. +double simulatePermeability(const Array &labels, + double voxel_length, double tau, + std::array force, int max_steps, + int interval, double tolerance, + const Utilities::MPI &comm); + +} // namespace rev_selection + +#endif diff --git a/tests/lbpm_permeability_simulator.cpp b/tests/lbpm_permeability_simulator.cpp index ad967104..e83818ad 100644 --- a/tests/lbpm_permeability_simulator.cpp +++ b/tests/lbpm_permeability_simulator.cpp @@ -11,6 +11,7 @@ #include "analysis/TwoPhase.h" #include "common/MPI.h" #include "models/MRTModel.h" +#include "analysis/rev_analysis.h" //#define WRITE_SURFACES /* @@ -41,14 +42,19 @@ int main(int argc, char **argv) comm.barrier(); ScaLBL_MRTModel MRT(rank,nprocs,comm); + ScaLBL_REV REV(rank,nprocs,comm); auto filename = argv[1]; MRT.ReadParams(filename); + REV.ReadParams(filename); MRT.SetDomain(); // this reads in the domain MRT.ReadInput(); MRT.Create(); // creating the model will create data structure to match the pore structure and allocate variables MRT.Initialize(); // initializing the model will set initial conditions for variables MRT.Run(); MRT.VelocityField(); + REV.DetRevAnalysis(MRT); + REV.StatRevAnalysis(MRT); + REV.sub_domain_select(MRT); cout << flush; } Utilities::shutdown(); diff --git a/tests/lbpm_rev_analysis.cpp b/tests/lbpm_rev_analysis.cpp new file mode 100644 index 00000000..ee6d5dae --- /dev/null +++ b/tests/lbpm_rev_analysis.cpp @@ -0,0 +1,96 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +#include "common/ScaLBL.h" +#include "common/Communication.h" +#include "common/Database.h" +#include "analysis/TwoPhase.h" +#include "common/MPI.h" +#include "analysis/rev_analysis.h" + +//#define WRITE_SURFACES + +/* + * Standalone representative elementary volume analysis + * James E. McClure 2013-2014 + */ + +using namespace std; + +namespace { +std::string missingHydraulicInputs(const std::string &filename) { + Database input_db(filename); + if (!input_db.keyExists("REV")) + return {}; + + auto rev_db = input_db.getDatabase("REV"); + if (!rev_db->getWithDefault("Hydraulic", false)) + return {}; + + std::string missing; + for (const char *key : {"SimulationParams", "VelocityPath", "VelocityTimestep"}) { + auto value = rev_db->getWithDefault(key, ""); + if (value.find_first_not_of(" \t\r\n") == std::string::npos) { + if (!missing.empty()) + missing += ", "; + missing += std::string("REV.") + key; + } + } + return missing; +} +} + + +int main(int argc, char **argv) +{ + if (argc < 2) { + std::cerr << "Usage: lbpm_rev_analysis params.db\n"; + return 1; + } + // Initialize MPI + Utilities::startup( argc, argv, false ); + Utilities::MPI comm( MPI_COMM_WORLD ); + int rank = comm.getRank(); + int nprocs = comm.getSize(); + int status = 0; + { + if (rank == 0){ + printf("********************************************************\n"); + printf("Running REV Analysis \n"); + printf("********************************************************\n"); + } + comm.barrier(); + + auto filename = argv[1]; + auto missing = missingHydraulicInputs(filename); + if (!missing.empty()) { + if (rank == 0) { + std::cerr << "lbpm_rev_analysis: Hydraulic = True requires you to provide " + << "SimulationParams, VelocityPath, and VelocityTimestep in the REV block.\n" + << "Missing or empty settings: " << missing + << ". Analysis was not run.\n"; + } + status = 1; + } else try { + ScaLBL_REV REV(rank,nprocs,comm); + REV.ReadParams(filename, false); + REV.SetDomain(); // this reads in the domain + REV.DetRevAnalysis(); + REV.StatRevAnalysis(); + REV.sub_domain_select(); + } catch (const std::exception &error) { + if (rank == 0) + std::cerr << "lbpm_rev_analysis: " << error.what() << '\n'; + status = 1; + } + cout << flush; + } + Utilities::shutdown(); + return status; +}