diff --git a/EXEs/cuda_mesh_state_report.cpp b/EXEs/cuda_mesh_state_report.cpp index 27d2aff..b26205a 100644 --- a/EXEs/cuda_mesh_state_report.cpp +++ b/EXEs/cuda_mesh_state_report.cpp @@ -1,5 +1,6 @@ #include "cuda/Cuda_mesh_state.hpp" #include "cuda/detail/Cuda_regular_geometry_cpu.hpp" +#include "cuda/detail/Cuda_regular_membrane_cpu.hpp" #include #include @@ -60,6 +61,35 @@ struct CaseResult } }; +struct MembraneCaseResult +{ + std::string name; + bool created = false; + bool cpuParity = true; + bool repeatable = true; + bool structuredDegeneracy = true; + bool recoverable = true; + bool permutationEqual = true; + double maxAbsError = 0.0; + std::vector firstValues; + std::vector firstCanonicalForces; + std::uint64_t warmedAllocations = 0; + DeviceStateReport activeReport; + DeviceStateReport finalReport; + DeviceStateError createError; + + bool pass() const + { + return created && cpuParity && repeatable && structuredDegeneracy && + recoverable && permutationEqual && + activeReport.phase == TransactionPhase::IdleAccepted && + finalReport.phase == TransactionPhase::Closed && + !finalReport.cleanupPending && finalReport.cleanupError.ok() && + finalReport.residentBytes == 0 && + finalReport.successfulAllocations == finalReport.successfulFrees; + } +}; + RegularMeshPack make_pack() { RegularMeshPack pack; @@ -106,6 +136,13 @@ RegularMeshPack make_pack() pack.shapeWeights[base + 26] = 1.0; } pack.parameters.kCurv = 1.0; + pack.parameters.uSurf = 2.0; + pack.parameters.uVol = 1.5; + pack.parameters.area0 = 1.0; + pack.parameters.area = 1.2; + pack.parameters.vol0 = 1.0; + pack.parameters.vol = 1.3; + pack.evaluatedFaceSpontaneousCurvature = {0.2}; return pack; } @@ -173,6 +210,14 @@ std::vector make_fixtures() production.shapeWeights[base + 13] = vScale; production.shapeWeights[base + 24] = -wScale; production.shapeWeights[base + 26] = wScale; + production.shapeWeights[base + 3 * 12] = -0.20; + production.shapeWeights[base + 3 * 12 + 3] = 0.20; + production.shapeWeights[base + 4 * 12] = -0.15; + production.shapeWeights[base + 4 * 12 + 4] = 0.15; + production.shapeWeights[base + 5 * 12 + 1] = -0.10; + production.shapeWeights[base + 5 * 12 + 5] = 0.10; + production.shapeWeights[base + 6 * 12 + 1] = -0.10; + production.shapeWeights[base + 6 * 12 + 5] = 0.10; } fixtures.push_back({"production_cpu", production, production.acceptedCoordinates}); @@ -203,6 +248,168 @@ bool same_bytes(const std::vector &left, left.size() * sizeof(double)) == 0); } +void append_values(std::vector &destination, + const std::vector &source) +{ + destination.insert(destination.end(), source.begin(), source.end()); +} + +std::vector membrane_values(const MembraneCandidateResult &result) +{ + std::vector values; + append_values(values, result.faceAreas); + append_values(values, result.faceVolumes); + append_values(values, result.faceBendingEnergies); + append_values(values, result.faceMeanCurvatures); + append_values(values, result.faceNormals); + append_values(values, result.occurrenceForces); + append_values(values, result.sampleSurfaceMeasures); + append_values(values, result.sampleMeanCurvatures); + append_values(values, result.sampleNormals); + append_values(values, result.sampleBendingEnergies); + values.push_back(result.totalArea); + values.push_back(result.totalVolume); + return values; +} + +std::vector membrane_values( + const detail::RegularMembraneCpuResult &result) +{ + std::vector values; + append_values(values, result.faceAreas); + append_values(values, result.faceVolumes); + append_values(values, result.faceBendingEnergies); + append_values(values, result.faceMeanCurvatures); + append_values(values, result.faceNormals); + append_values(values, result.occurrenceForces); + append_values(values, result.sampleSurfaceMeasures); + append_values(values, result.sampleMeanCurvatures); + append_values(values, result.sampleNormals); + append_values(values, result.sampleBendingEnergies); + values.push_back(result.totalArea); + values.push_back(result.totalVolume); + return values; +} + +double max_abs_error(const std::vector &actual, + const std::vector &expected) +{ + if (actual.size() != expected.size()) + return std::numeric_limits::infinity(); + double error = 0.0; + for (std::size_t index = 0; index < actual.size(); ++index) + error = std::max(error, std::abs(actual[index] - expected[index])); + return error; +} + +std::vector canonical_occurrence_forces( + const RegularMeshPack &pack, + const std::vector &occurrenceForces) +{ + std::vector result( + static_cast(pack.vertexCount) * 9U, 0.0); + for (std::size_t evaluated = 0; + evaluated < static_cast(pack.evaluatedFaceCount); + ++evaluated) + for (std::size_t local = 0; local < kRegularControlCount; ++local) + { + const std::size_t occurrence = + evaluated * kRegularControlCount + local; + const std::size_t source = static_cast( + pack.oneRingSourceIds[occurrence]); + for (std::size_t component = 0; component < 9; ++component) + result[source * 9U + component] = + occurrenceForces[occurrence * 9U + component]; + } + return result; +} + +MembraneCaseResult run_membrane_case(const GeometryFixture &fixture, + int device, int iterations) +{ + MembraneCaseResult result; + result.name = fixture.name; + DeviceStateConfig config; + config.deviceOrdinal = device; + auto created = create_cuda_mesh_state(fixture.pack, config); + result.createError = created.report.error; + if (!created.ok()) + { + result.finalReport = created.report; + return result; + } + result.created = true; + result.warmedAllocations = created.state->report().successfulAllocations; + const detail::RegularMembraneCpuResult expected = + detail::evaluate_regular_membrane_cpu(fixture.pack, + fixture.candidate); + const bool expectDegenerate = fixture.name == "degenerate"; + result.structuredDegeneracy = expectDegenerate ? !expected.ok() + : expected.ok(); + const std::vector expectedValues = membrane_values(expected); + for (int iteration = 0; iteration < iterations; ++iteration) + { + const std::uint64_t generation = + created.state->report().residentGenerations.acceptedCoordinates + + 100U + static_cast(iteration); + if (!created.state->prepare_candidate(fixture.candidate, + generation).ok()) + { + result.recoverable = false; + break; + } + const MembraneCandidateResult membrane = + created.state->compute_candidate_membrane(); + if (expectDegenerate) + { + result.structuredDegeneracy = result.structuredDegeneracy && + !membrane.ok() && + membrane.status == MembraneStatusCode::DegenerateSample && + membrane.error.code == DeviceStateErrorCode::CandidateFailed; + if (!created.state->recover().ok()) + { + result.recoverable = false; + break; + } + continue; + } + if (!membrane.ok()) + { + result.cpuParity = false; + result.recoverable = false; + break; + } + const std::vector actualValues = membrane_values(membrane); + result.maxAbsError = std::max( + result.maxAbsError, + max_abs_error(actualValues, expectedValues)); + if (iteration == 0) + { + result.firstValues = actualValues; + result.firstCanonicalForces = canonical_occurrence_forces( + fixture.pack, membrane.occurrenceForces); + } + else + result.repeatable = result.repeatable && + same_bytes(result.firstValues, actualValues); + if (!created.state->rollback().ok()) + { + result.recoverable = false; + break; + } + } + result.cpuParity = result.cpuParity && + (expectDegenerate || result.maxAbsError <= kTolerance); + result.activeReport = created.state->report(); + result.recoverable = result.recoverable && + result.activeReport.phase == TransactionPhase::IdleAccepted && + result.activeReport.successfulAllocations == result.warmedAllocations; + const DeviceStateError closed = created.state->close(); + result.finalReport = created.state->report(); + result.recoverable = result.recoverable && closed.ok(); + return result; +} + CaseResult run_case(const GeometryFixture &fixture, int device, int iterations) { CaseResult result; @@ -317,8 +524,9 @@ int main(int argc, char **argv) { const int device = integer_argument(argv, argc, "--device", 0); const int iterations = integer_argument(argv, argc, "--iterations", 20); + const std::vector fixtures = make_fixtures(); std::vector cases; - for (const GeometryFixture &fixture : make_fixtures()) + for (const GeometryFixture &fixture : fixtures) { cases.push_back(run_case(fixture, device, iterations)); if (!cases.back().created) @@ -335,6 +543,10 @@ int main(int argc, char **argv) return report.compiled ? 77 : 0; } } + std::vector membraneCases; + for (const GeometryFixture &fixture : fixtures) + membraneCases.push_back(run_membrane_case(fixture, device, + iterations)); const CaseResult &natural = cases[0]; CaseResult &permuted = cases[1]; @@ -345,6 +557,9 @@ int main(int argc, char **argv) sizeof(double)) == 0 && std::memcmp(&natural.firstTotalVolume, &permuted.firstTotalVolume, sizeof(double)) == 0; + membraneCases[1].permutationEqual = + same_bytes(membraneCases[0].firstCanonicalForces, + membraneCases[1].firstCanonicalForces); bool pass = true; bool geometryRepeatable = true; @@ -354,8 +569,12 @@ int main(int argc, char **argv) bool closed = true; bool cleanupPending = false; double geometryMaxAbsError = 0.0; + double membraneMaxAbsError = 0.0; + bool membraneRepeatable = true; + bool membraneDegeneracyHandled = true; std::uint64_t commits = 0, rollbacks = 0, allocationEpoch = 0; std::uint64_t transactionEpoch = 0, warmedAllocations = 0; + std::uint64_t membraneTransactions = 0; std::uint64_t finalAllocations = 0, successfulFrees = 0; std::size_t residentBytes = 0, finalResidentBytes = 0; std::size_t memoryBudgetBytes = 0; @@ -394,7 +613,47 @@ int main(int argc, char **argv) transfers[reason].completedBytes += source.completedBytes; } } + for (const MembraneCaseResult &item : membraneCases) + { + pass = pass && item.pass(); + membraneRepeatable = membraneRepeatable && item.repeatable; + membraneDegeneracyHandled = membraneDegeneracyHandled && + item.structuredDegeneracy; + membraneMaxAbsError = std::max(membraneMaxAbsError, + item.maxAbsError); + membraneTransactions += item.activeReport.transactionEpoch; + transactionEpoch += item.activeReport.transactionEpoch; + allocationEpoch += item.activeReport.allocationEpoch; + warmedAllocations += item.warmedAllocations; + finalAllocations += item.finalReport.successfulAllocations; + successfulFrees += item.finalReport.successfulFrees; + residentBytes = std::max(residentBytes, + item.activeReport.residentBytes); + finalResidentBytes += item.finalReport.residentBytes; + noWarmAllocations = noWarmAllocations && + item.activeReport.successfulAllocations == item.warmedAllocations; + allocationFreeBalance = allocationFreeBalance && + item.finalReport.successfulAllocations == + item.finalReport.successfulFrees; + closed = closed && + item.finalReport.phase == TransactionPhase::Closed; + cleanupPending = cleanupPending || + item.finalReport.cleanupPending; + for (std::size_t reason = 0; reason < transfers.size(); ++reason) + { + const TransferCounter &source = + item.activeReport.transfers[reason]; + transfersComplete = transfersComplete && + source.attemptedOperations == source.completedOperations && + source.attemptedBytes == source.completedBytes; + transfers[reason].completedOperations += + source.completedOperations; + transfers[reason].completedBytes += source.completedBytes; + } + } pass = pass && geometryMaxAbsError <= kTolerance && geometryRepeatable && + membraneMaxAbsError <= kTolerance && membraneRepeatable && + membraneDegeneracyHandled && transfersComplete && noWarmAllocations && allocationFreeBalance && closed && !cleanupPending && finalResidentBytes == 0; @@ -403,7 +662,9 @@ int main(int argc, char **argv) << ",\"device_ordinal\":" << device << ",\"iterations\":" << iterations << ",\"case_count\":" << cases.size() - << ",\"total_transactions\":" << commits + rollbacks + << ",\"total_transactions\":" << transactionEpoch + << ",\"membrane_transactions\":" + << membraneTransactions << ",\"commits\":" << commits << ",\"rollbacks\":" << rollbacks << ",\"allocation_epoch\":" << allocationEpoch @@ -419,6 +680,12 @@ int main(int argc, char **argv) << geometryMaxAbsError << ",\"geometry_repeatable\":" << (geometryRepeatable ? "true" : "false") + << ",\"membrane_max_abs_error\":" + << membraneMaxAbsError + << ",\"membrane_repeatable\":" + << (membraneRepeatable ? "true" : "false") + << ",\"membrane_degeneracy_handled\":" + << (membraneDegeneracyHandled ? "true" : "false") << ",\"resident_bytes\":" << residentBytes << ",\"final_resident_bytes\":" << finalResidentBytes << ",\"closed\":" << (closed ? "true" : "false") @@ -448,6 +715,26 @@ int main(int argc, char **argv) << ",\"permutation_equal\":" << (item.permutationEqual ? "true" : "false") << '}'; } + std::cout << "},\"membrane_cases\":{"; + for (std::size_t index = 0; index < membraneCases.size(); ++index) + { + if (index) + std::cout << ','; + const MembraneCaseResult &item = membraneCases[index]; + std::cout << json_string(item.name) + << ":{\"pass\":" << (item.pass() ? "true" : "false") + << ",\"cpu_parity\":" + << (item.cpuParity ? "true" : "false") + << ",\"repeatable\":" + << (item.repeatable ? "true" : "false") + << ",\"structured_degeneracy\":" + << (item.structuredDegeneracy ? "true" : "false") + << ",\"recoverable\":" + << (item.recoverable ? "true" : "false") + << ",\"permutation_equal\":" + << (item.permutationEqual ? "true" : "false") + << ",\"max_abs_error\":" << item.maxAbsError << '}'; + } std::cout << "},\"transfers\":{"; for (std::size_t reason = 0; reason < transfers.size(); ++reason) { diff --git a/Makefile b/Makefile index d57f602..c8cae63 100644 --- a/Makefile +++ b/Makefile @@ -355,7 +355,9 @@ cuda_backend_stub_report: include/cuda/Cuda_backend.hpp \ cuda_mesh_state_report: include/cuda/Cuda_mesh_state.hpp \ include/cuda/detail/Cuda_mesh_state_core.hpp \ include/cuda/detail/Cuda_regular_geometry_cpu.hpp \ + include/cuda/detail/Cuda_regular_membrane_cpu.hpp \ src/cuda/Cuda_regular_geometry_cpu.cpp \ + src/cuda/Cuda_regular_membrane_cpu.cpp \ src/cuda/Cuda_mesh_state_common.cpp src/cuda/Cuda_mesh_state.cu \ EXEs/cuda_mesh_state_report.cpp @command -v $(CUDA_NVCC) >/dev/null 2>&1 || \ @@ -364,6 +366,7 @@ cuda_mesh_state_report: include/cuda/Cuda_mesh_state.hpp \ -arch=$(CUDA_COMPUTE_ARCH) -code=$(CUDA_SM_CODE) \ -ccbin=$(CUDA_HOST_CXX) -Iinclude \ src/cuda/Cuda_regular_geometry_cpu.cpp \ + src/cuda/Cuda_regular_membrane_cpu.cpp \ src/cuda/Cuda_mesh_state_common.cpp src/cuda/Cuda_mesh_state.cu \ EXEs/cuda_mesh_state_report.cpp \ -lcudart -o $(CUDA_MESH_STATE_REPORT) @@ -371,11 +374,14 @@ cuda_mesh_state_report: include/cuda/Cuda_mesh_state.hpp \ cuda_mesh_state_stub_report: include/cuda/Cuda_mesh_state.hpp \ include/cuda/detail/Cuda_regular_geometry_cpu.hpp \ + include/cuda/detail/Cuda_regular_membrane_cpu.hpp \ src/cuda/Cuda_regular_geometry_cpu.cpp \ + src/cuda/Cuda_regular_membrane_cpu.cpp \ src/cuda/Cuda_mesh_state_common.cpp src/cuda/Cuda_mesh_state_stub.cpp \ EXEs/cuda_mesh_state_report.cpp $(CXX) -std=$(CXX_STD) -O3 -Iinclude \ src/cuda/Cuda_regular_geometry_cpu.cpp \ + src/cuda/Cuda_regular_membrane_cpu.cpp \ src/cuda/Cuda_mesh_state_common.cpp src/cuda/Cuda_mesh_state_stub.cpp \ EXEs/cuda_mesh_state_report.cpp -o $(CUDA_MESH_STATE_STUB_REPORT) @echo "Finished non-CUDA mesh-state stub report build, $(CUDA_MESH_STATE_STUB_REPORT)." diff --git a/analysis/cuda_mesh_state_report_rtx4050.json b/analysis/cuda_mesh_state_report_rtx4050.json index f7542e2..a31ae15 100644 --- a/analysis/cuda_mesh_state_report_rtx4050.json +++ b/analysis/cuda_mesh_state_report_rtx4050.json @@ -1,5 +1,5 @@ { - "allocation_epoch": 6, + "allocation_epoch": 12, "allocation_free_balance": true, "available": true, "build": { @@ -27,7 +27,7 @@ "compiled": true, "cuda_required": true, "device_ordinal": 0, - "final_allocations": 138, + "final_allocations": 384, "final_resident_bytes": 0, "geometry_cases": { "boundary_ghost": { @@ -94,24 +94,116 @@ "name": "NVIDIA GeForce RTX 4050 Laptop GPU" }, "iterations": 20, + "membrane_cases": { + "boundary_ghost": { + "cpu_parity": true, + "max_abs_error": 0, + "pass": true, + "permutation_equal": true, + "recoverable": true, + "repeatable": true, + "structured_degeneracy": true + }, + "curved": { + "cpu_parity": true, + "max_abs_error": 0, + "pass": true, + "permutation_equal": true, + "recoverable": true, + "repeatable": true, + "structured_degeneracy": true + }, + "degenerate": { + "cpu_parity": true, + "max_abs_error": 0, + "pass": true, + "permutation_equal": true, + "recoverable": true, + "repeatable": true, + "structured_degeneracy": true + }, + "natural": { + "cpu_parity": true, + "max_abs_error": 0, + "pass": true, + "permutation_equal": true, + "recoverable": true, + "repeatable": true, + "structured_degeneracy": true + }, + "permuted": { + "cpu_parity": true, + "max_abs_error": 0, + "pass": true, + "permutation_equal": true, + "recoverable": true, + "repeatable": true, + "structured_degeneracy": true + }, + "production_cpu": { + "cpu_parity": true, + "max_abs_error": 3.1086244689504383e-15, + "pass": true, + "permutation_equal": true, + "recoverable": true, + "repeatable": true, + "structured_degeneracy": true + } + }, + "membrane_degeneracy_handled": true, + "membrane_max_abs_error": 3.1086244689504383e-15, + "membrane_repeatable": true, + "membrane_transactions": 120, "memory_budget_bytes": 2659713024, "no_warm_allocations": true, - "resident_bytes": 5504, + "resident_bytes": 7104, "rollbacks": 60, "status": "pass", - "successful_frees": 138, - "total_transactions": 120, - "transaction_epoch": 120, + "successful_frees": 384, + "total_transactions": 240, + "transaction_epoch": 240, "transfers": { - "accepted_coordinates": {"bytes": 3456, "operations": 12}, - "candidate_coordinates": {"bytes": 34560, "operations": 120}, - "candidate_geometry": {"bytes": 0, "operations": 0}, - "geometry_diagnostics": {"bytes": 4640, "operations": 480}, - "numerical_plan": {"bytes": 12672, "operations": 18}, - "parameters": {"bytes": 822, "operations": 18}, - "reference_coordinates": {"bytes": 1728, "operations": 6}, - "topology": {"bytes": 1742, "operations": 54} + "accepted_coordinates": { + "bytes": 6912, + "operations": 24 + }, + "candidate_coordinates": { + "bytes": 69120, + "operations": 240 + }, + "candidate_geometry": { + "bytes": 0, + "operations": 0 + }, + "candidate_membrane": { + "bytes": 0, + "operations": 0 + }, + "geometry_diagnostics": { + "bytes": 4640, + "operations": 480 + }, + "membrane_diagnostics": { + "bytes": 132160, + "operations": 1440 + }, + "numerical_plan": { + "bytes": 25344, + "operations": 36 + }, + "parameters": { + "bytes": 1644, + "operations": 36 + }, + "reference_coordinates": { + "bytes": 3456, + "operations": 12 + }, + "topology": { + "bytes": 3484, + "operations": 108 + } }, "transfers_complete": true, - "warmed_allocations": 138 + "warmed_allocations": 384 } diff --git a/docs/cuda_end_to_end_residency_force_scatter_implementation_plan.md b/docs/cuda_end_to_end_residency_force_scatter_implementation_plan.md index 237a02a..f13f370 100644 --- a/docs/cuda_end_to_end_residency_force_scatter_implementation_plan.md +++ b/docs/cuda_end_to_end_residency_force_scatter_implementation_plan.md @@ -740,6 +740,21 @@ remains subject to the required pull-request review and owner merge approval. Exit evidence: per-sample, per-face, per-occurrence, energy, curvature, and normal parity against the actual CPU formula. +Implementation status (2026-08-03): +`codex/cuda-step5-regular-membrane-formula` implements the complete regular +bending, global-area-constraint, and global-volume-constraint sample algebra +as a candidate-only device operation. It consumes all seven packed weighted +rows and emits sample surface measures/curvatures/normals/bending energies, +face geometry/energy/curvature/unit normals, and canonical 12x9 +per-occurrence force contributions. An independent packed CPU oracle is bound +directly to `Mesh::element_energy_force_regular()` for every face observable +and force component on a curved production mesh. The native RTX report covers +the six frozen fixtures for 20 repeats, including source-order permutation and +structured degenerate-sample recovery, without warmed allocations. Scatter, +vertex-force writes, `Mesh` publication, and production routing remain outside +this slice. Final status remains subject to pull-request review and owner merge +approval. + ### Step 6 / PR 6: Deterministic source-keyed scatter - Implement the reviewed incidence-based, fixed-order reduction into nine diff --git a/docs/cuda_mesh_state.md b/docs/cuda_mesh_state.md index e10a652..d7f26b1 100644 --- a/docs/cuda_mesh_state.md +++ b/docs/cuda_mesh_state.md @@ -1,10 +1,11 @@ -# Persistent CUDA Mesh State, Transactions, And Geometry +# Persistent CUDA Mesh State, Geometry, And Regular Membrane Candidates -This document covers Steps 3 and 4 of the end-to-end CUDA residency program. +This document covers Steps 3 through 5 of the end-to-end CUDA residency program. Step 3 added persistent storage and transaction control. Step 4 adds only the eligible regular-face area/volume calculation and deterministic global -area/volume reductions. There is still no membrane-force formula, force -scatter, evaluator routing, or host `Mesh` publication. +area/volume reductions. Step 5 adds the complete regular membrane formula into +candidate buffers. There is still no force scatter, evaluator routing, or host +`Mesh` publication. ## Authority and ownership @@ -14,8 +15,8 @@ translation unit; ordinary builds discover the structured non-CUDA stub. Accepted host state remains authoritative until later routing steps. The persistent groups are topology/incidence, numerical plan, parameters, -accepted/previous/candidate coordinates, reference coordinates, and candidate -geometry/status buffers. Every +accepted/previous/candidate coordinates, reference coordinates, candidate +geometry/status buffers, and membrane diagnostic/output buffers. Every group is keyed by the corresponding `MeshPackGenerations` value. Equal generations cause no allocation, copy, or synchronization. A changed group is allocated geometrically, copied into staging, synchronized, and installed only @@ -66,6 +67,26 @@ The second expression intentionally retains the existing legacy first- component `dot_row` behavior. A device status rejects invalid indices or nonfinite results before the candidate becomes validated. +`compute_candidate_membrane()` is an alternative candidate operation from +`CandidatePrepared`. It consumes all seven weighted rows and the packed +per-face spontaneous curvature plus global bending/area/volume parameters. It +emits, without scatter: + +- per-sample surface measure, mean curvature, unit normal, and bending energy; +- per-face area, legacy volume, bending energy, integrated mean curvature, and + normalized face normal; and +- 12 occurrences x 9 components ordered as bending xyz, area xyz, volume xyz. + +The device formula follows `Mesh::element_energy_force_regular()` including +the current half-quadrature convention, dual-basis derivatives, bending +gradient matrices, and global constraint factors. A structured status records +the first invalid source, degenerate sample, nonfinite intermediate, or +nonfinite output together with evaluated-face and sample indices. Such a +candidate enters `Failed`; `recover()` returns to the unchanged accepted state. +Unlike Step 4 geometry, a zero surface measure is invalid for Step 5 because +the force formula requires its inverse. No source incidence plan is consumed +and no vertex-indexed force buffer is written. + Commit rotates candidate into accepted, accepted into previous, and the old previous slot into reusable candidate storage while advancing the accepted coordinate generation without a mesh-sized copy. Rollback from any live @@ -78,10 +99,10 @@ searches; it is not hidden or misclassified here. Each transfer records attempted/completed operation and byte counts under a stable reason: topology, numerical plan, parameters, accepted coordinates, -reference coordinates, candidate coordinates, candidate geometry storage, or -geometry diagnostics. Step 4 copies face/global outputs to the host only for -the explicit comparison API and report; those mesh-sized diagnostic copies are -classified and are not a production route. Allocation and +reference coordinates, candidate coordinates, candidate geometry or membrane +storage, and geometry or membrane diagnostics. Steps 4 and 5 copy outputs to +the host only for their explicit comparison APIs and reports; those diagnostic +copies are classified and are not a production route. Allocation and transaction epochs, live bytes, allocation/free counts, synchronizations, roles, generations, and the last outcome are also reportable. `lastDirtyGroups` is reset at the start of each candidate-upload decision. A @@ -107,7 +128,11 @@ curved production regular mesh at the `1.0e-12` gate. Injected kernel, diagnostic-copy, and synchronization failures and injected nonzero status, nonfinite face output, negative area, or nonfinite totals remain recoverable without validating or changing accepted -state. The focused state suite passes 27/27 tests. +state. Step 5 additionally rejects injected kernel, diagnostic-copy, +synchronization, status, nonfinite-force, and zero-surface failures while +preserving accepted bytes and recoverability. Its independent packed CPU +oracle is compared directly with the production regular formula at the +per-face and per-occurrence levels. The explicit native proof builds and runs with: @@ -118,16 +143,21 @@ python3 scripts/run_cuda_mesh_state_report.py \ The native report executes the natural, permuted, sample-varying curved, boundary/ghost, degenerate, and production-formula CPU-oracle fixtures on the -actual CUDA kernel. Each case runs 20 candidate geometry transactions and has -separate parity, structural, and bitwise-repeatability fields; the runner -rejects a report with any missing or false case field. The aggregate requires -`geometry_max_abs_error <= 1.0e-12` plus bitwise repeatability. It also closes +actual CUDA kernels. Each case runs 20 geometry and 20 membrane candidate +transactions. Membrane evidence binds every sample diagnostic, face +observable, and 12x9 occurrence contribution; the degenerate case must report +the structured failure and recover on every repeat. The runner rejects a +report with any missing or false case field. The aggregate requires both +geometry and membrane maximum error at or below `1.0e-12` plus bitwise +repeatability. It also closes explicitly before declaring success and requires `Closed`, zero cleanup debt, zero final resident bytes, and exact allocation/free balance. Geometry parity, repeatability, and teardown are therefore part of the RTX pass predicate, not informational fields or destructor side effects after reporting. The recorded -RTX proof covers 120 transactions with maximum error -`1.3877787807814457e-17`, no warmed allocations, and 138/138 balanced frees. +RTX proof covers 240 transactions. Its membrane maximum absolute error is +`3.1086244689504383e-15`; geometry remains +`1.3877787807814457e-17`. It has no warmed allocations and exact allocation/ +free balance. The non-CUDA contract is independently runnable with `--stub`. Native evidence for the development RTX machine is recorded in diff --git a/include/cuda/Cuda_mesh_state.hpp b/include/cuda/Cuda_mesh_state.hpp index 8c09ed5..35588c8 100644 --- a/include/cuda/Cuda_mesh_state.hpp +++ b/include/cuda/Cuda_mesh_state.hpp @@ -78,6 +78,8 @@ enum class TransferReason : std::size_t CandidateCoordinates, CandidateGeometry, GeometryDiagnostics, + CandidateMembrane, + MembraneDiagnostics, Count, }; @@ -142,6 +144,41 @@ struct GeometryCandidateResult bool ok() const noexcept { return error.ok(); } }; +enum class MembraneStatusCode : std::int32_t +{ + None = 0, + InvalidSource = 1, + DegenerateSample = 2, + NonFiniteIntermediate = 3, + NonFiniteOutput = 4, +}; + +struct MembraneCandidateResult +{ + std::vector faceAreas; + std::vector faceVolumes; + std::vector faceBendingEnergies; + std::vector faceMeanCurvatures; + std::vector faceNormals; + std::vector occurrenceForces; + std::vector sampleSurfaceMeasures; + std::vector sampleMeanCurvatures; + std::vector sampleNormals; + std::vector sampleBendingEnergies; + double totalArea = 0.0; + double totalVolume = 0.0; + std::uint64_t coordinateGeneration = 0; + MembraneStatusCode status = MembraneStatusCode::None; + std::uint64_t failedEvaluatedFace = 0; + std::uint32_t failedSample = 0; + DeviceStateError error; + + bool ok() const noexcept + { + return error.ok() && status == MembraneStatusCode::None; + } +}; + struct CudaMeshStateResult; class CudaMeshState final @@ -158,6 +195,7 @@ class CudaMeshState final const std::vector &coordinates, std::uint64_t generation); GeometryCandidateResult compute_candidate_geometry(); + MembraneCandidateResult compute_candidate_membrane(); DeviceStateError mark_computing(); DeviceStateError mark_validated(); DeviceStateError commit(); diff --git a/include/cuda/detail/Cuda_mesh_state_core.hpp b/include/cuda/detail/Cuda_mesh_state_core.hpp index 0406460..fc26c4d 100644 --- a/include/cuda/detail/Cuda_mesh_state_core.hpp +++ b/include/cuda/detail/Cuda_mesh_state_core.hpp @@ -37,6 +37,32 @@ struct GeometryLaunch std::uint64_t evaluatedFaceCount = 0; }; +struct MembraneLaunch +{ + DeviceBufferHandle evaluatedFaceIds = 0; + DeviceBufferHandle oneRingSourceIds = 0; + DeviceBufferHandle spontaneousCurvatures = 0; + DeviceBufferHandle packedParameters = 0; + DeviceBufferHandle quadratureCoefficients = 0; + DeviceBufferHandle shapeWeights = 0; + DeviceBufferHandle coordinates = 0; + DeviceBufferHandle faceAreas = 0; + DeviceBufferHandle faceVolumes = 0; + DeviceBufferHandle geometryTotals = 0; + DeviceBufferHandle faceBendingEnergies = 0; + DeviceBufferHandle faceMeanCurvatures = 0; + DeviceBufferHandle faceNormals = 0; + DeviceBufferHandle occurrenceForces = 0; + DeviceBufferHandle sampleSurfaceMeasures = 0; + DeviceBufferHandle sampleMeanCurvatures = 0; + DeviceBufferHandle sampleNormals = 0; + DeviceBufferHandle sampleBendingEnergies = 0; + DeviceBufferHandle statusDiagnostics = 0; + std::uint64_t vertexCount = 0; + std::uint64_t faceCount = 0; + std::uint64_t evaluatedFaceCount = 0; +}; + DriverStatus release_retryable_handle( DeviceBufferHandle &handle, const std::function &release); @@ -66,6 +92,7 @@ struct DeviceOperations std::function copyDeviceToHost; std::function computeGeometry; + std::function computeMembrane; std::function synchronize; }; @@ -82,6 +109,7 @@ class MeshStateCore final DeviceStateError prepare_candidate(const std::vector &coordinates, std::uint64_t generation); GeometryCandidateResult compute_candidate_geometry(); + MembraneCandidateResult compute_candidate_membrane(); DeviceStateError mark_computing(); DeviceStateError mark_validated(); DeviceStateError commit(); diff --git a/include/cuda/detail/Cuda_regular_membrane_cpu.hpp b/include/cuda/detail/Cuda_regular_membrane_cpu.hpp new file mode 100644 index 0000000..43356dd --- /dev/null +++ b/include/cuda/detail/Cuda_regular_membrane_cpu.hpp @@ -0,0 +1,59 @@ +#ifndef SLIMED_CUDA_REGULAR_MEMBRANE_CPU_HPP +#define SLIMED_CUDA_REGULAR_MEMBRANE_CPU_HPP + +#include "cuda/Cuda_mesh_pack.hpp" + +#include +#include +#include + +namespace slimed::cuda_residency::detail +{ + +enum class RegularMembraneStatus : std::int32_t +{ + None = 0, + InvalidSource = 1, + DegenerateSample = 2, + NonFiniteIntermediate = 3, + NonFiniteOutput = 4, +}; + +const char *regular_membrane_status_name(RegularMembraneStatus status) noexcept; + +/** + * Host-side oracle for the packed regular membrane contract. + * + * Face arrays use declared face IDs. Sample arrays and occurrence forces use + * evaluated-face order. Occurrence force components are ordered as + * [bending xyz, area xyz, volume xyz]. No source-keyed scatter is performed. + */ +struct RegularMembraneCpuResult +{ + std::vector faceAreas; + std::vector faceVolumes; + std::vector faceBendingEnergies; + std::vector faceMeanCurvatures; + std::vector faceNormals; + std::vector occurrenceForces; + std::vector sampleSurfaceMeasures; + std::vector sampleMeanCurvatures; + std::vector sampleNormals; + std::vector sampleBendingEnergies; + double totalArea = 0.0; + double totalVolume = 0.0; + RegularMembraneStatus status = RegularMembraneStatus::None; + std::uint64_t failedEvaluatedFace = 0; + std::uint32_t failedSample = 0; + std::string message; + + bool ok() const noexcept { return status == RegularMembraneStatus::None; } +}; + +RegularMembraneCpuResult evaluate_regular_membrane_cpu( + const RegularMeshPack &pack, + const std::vector &coordinates); + +} // namespace slimed::cuda_residency::detail + +#endif diff --git a/scripts/run_cuda_mesh_state_report.py b/scripts/run_cuda_mesh_state_report.py index 7adcc99..088a5d7 100644 --- a/scripts/run_cuda_mesh_state_report.py +++ b/scripts/run_cuda_mesh_state_report.py @@ -22,6 +22,7 @@ "degenerate", "production_cpu", } +REQUIRED_MEMBRANE_CASES = REQUIRED_GEOMETRY_CASES def repo_root() -> Path: @@ -128,6 +129,35 @@ def geometry_complete(report: dict[str, object]) -> bool: ) +def membrane_complete(report: dict[str, object]) -> bool: + error = report.get("membrane_max_abs_error") + cases = report.get("membrane_cases") + if not isinstance(cases, dict) or set(cases) != REQUIRED_MEMBRANE_CASES: + return False + for name in REQUIRED_MEMBRANE_CASES: + case = cases.get(name) + if not isinstance(case, dict): + return False + case_error = case.get("max_abs_error") + if not ( + case.get("pass") is True + and case.get("cpu_parity") is True + and case.get("repeatable") is True + and case.get("structured_degeneracy") is True + and case.get("recoverable") is True + and case.get("permutation_equal") is True + and isinstance(case_error, (int, float)) + and case_error <= 1.0e-12 + ): + return False + return ( + report.get("membrane_repeatable") is True + and report.get("membrane_degeneracy_handled") is True + and isinstance(error, (int, float)) + and error <= 1.0e-12 + ) + + def run(args: argparse.Namespace) -> int: root = Path(args.root).resolve() if args.root else repo_root() make = executable(args.make, "make") @@ -175,6 +205,11 @@ def run(args: argparse.Namespace) -> int: "mesh-state report claimed pass without geometry parity and repeatability\n" ) return 1 + if not membrane_complete(report): + sys.stderr.write( + "mesh-state report claimed pass without complete membrane parity, repeatability, and recovery\n" + ) + return 1 if not teardown_complete(report): sys.stderr.write( "mesh-state report claimed pass without complete, balanced teardown\n" @@ -193,7 +228,13 @@ def run(args: argparse.Namespace) -> int: } if not args.stub: report["gpu"] = gpu_inventory() - print(json.dumps(report, indent=2, sort_keys=True)) + encoded = json.dumps(report, indent=2, sort_keys=True) + if args.output: + output = Path(args.output) + if not output.is_absolute(): + output = root / output + output.write_text(encoded + "\n") + print(encoded) if execution.returncode == NO_CUDA_EXIT_CODE: return NO_CUDA_EXIT_CODE if args.require_cuda and not args.stub else 0 return execution.returncode @@ -202,6 +243,7 @@ def run(args: argparse.Namespace) -> int: def parser() -> argparse.ArgumentParser: result = argparse.ArgumentParser(description=__doc__) result.add_argument("--root") + result.add_argument("--output") result.add_argument("--make") result.add_argument("--nvcc") result.add_argument("--host-cxx") diff --git a/src/cuda/Cuda_mesh_state.cu b/src/cuda/Cuda_mesh_state.cu index bf48899..4a60bf9 100644 --- a/src/cuda/Cuda_mesh_state.cu +++ b/src/cuda/Cuda_mesh_state.cu @@ -14,6 +14,75 @@ namespace constexpr double kLegacyVolumeQuadratureFactor = 0.16666666666; +__device__ double dot3(const double left[3], const double right[3]) +{ + return left[0] * right[0] + left[1] * right[1] + + left[2] * right[2]; +} + +__device__ void cross3(const double left[3], const double right[3], + double result[3]) +{ + result[0] = left[1] * right[2] - left[2] * right[1]; + result[1] = left[2] * right[0] - left[0] * right[2]; + result[2] = left[0] * right[1] - left[1] * right[0]; +} + +__device__ void add3(const double left[3], const double right[3], + double result[3]) +{ + for (int axis = 0; axis < 3; ++axis) + result[axis] = left[axis] + right[axis]; +} + +__device__ void linear3(const double left[3], double leftFactor, + const double right[3], double rightFactor, + double result[3]) +{ + for (int axis = 0; axis < 3; ++axis) + result[axis] = leftFactor * left[axis] + + rightFactor * right[axis]; +} + +__device__ bool finite3(const double value[3]) +{ + return isfinite(value[0]) && isfinite(value[1]) && isfinite(value[2]); +} + +__device__ void add_outer_scaled(double matrix[3][3], + const double left[3], + const double right[3], double factor) +{ + for (int row = 0; row < 3; ++row) + for (int column = 0; column < 3; ++column) + matrix[row][column] += + factor * left[row] * right[column]; +} + +__device__ void transpose_multiply3(const double matrix[3][3], + const double value[3], + double result[3]) +{ + for (int column = 0; column < 3; ++column) + { + result[column] = 0.0; + for (int row = 0; row < 3; ++row) + result[column] += value[row] * matrix[row][column]; + } +} + +__device__ void record_membrane_status(std::int32_t *diagnostics, + std::int32_t code, + std::uint64_t evaluated, + std::uint32_t sample) +{ + if (atomicCAS(diagnostics, 0, code) == 0) + { + diagnostics[1] = static_cast(evaluated); + diagnostics[2] = static_cast(sample); + } +} + __global__ void regular_geometry_kernel( const std::int32_t *evaluatedFaceIds, const std::int32_t *oneRingSourceIds, @@ -106,6 +175,335 @@ __global__ void deterministic_geometry_reduction_kernel( totals[1] = volume; } +__global__ void regular_membrane_kernel( + const std::int32_t *evaluatedFaceIds, + const std::int32_t *oneRingSourceIds, + const double *spontaneousCurvatures, + const PackedRegularParameters *parameters, + const double *quadratureCoefficients, + const double *shapeWeights, + const double *coordinates, + double *faceAreas, + double *faceVolumes, + double *faceBendingEnergies, + double *faceMeanCurvatures, + double *faceNormals, + double *occurrenceForces, + double *sampleSurfaceMeasures, + double *sampleMeanCurvatures, + double *sampleNormals, + double *sampleBendingEnergies, + std::int32_t *diagnostics, + std::uint64_t vertexCount, + std::uint64_t faceCount, + std::uint64_t evaluatedFaceCount) +{ + const std::uint64_t evaluated = + static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + if (evaluated >= evaluatedFaceCount) + return; + const std::int32_t faceId = evaluatedFaceIds[evaluated]; + if (faceId < 0 || static_cast(faceId) >= faceCount) + { + record_membrane_status(diagnostics, 1, evaluated, 0); + return; + } + + const PackedRegularParameters parameter = *parameters; + const double uSurfPerArea = + parameter.uSurf == 0.0 || parameter.area0 == 0.0 + ? 0.0 + : parameter.uSurf / parameter.area0; + const double uVol = + parameter.uVol == 0.0 || parameter.vol0 == 0.0 + ? 0.0 + : parameter.uVol / parameter.vol0; + const double areaFactor = + uSurfPerArea * (parameter.area - parameter.area0); + const double volumeFactor = + uVol * (parameter.vol - parameter.vol0) / 3.0; + double area = 0.0; + double volume = 0.0; + double faceBending = 0.0; + double faceMean = 0.0; + double accumulatedNormal[3]{}; + + for (std::uint32_t sample = 0; sample < kQuadratureSampleCount; ++sample) + { + double rows[kShapeRowCount][3]{}; + for (std::uint32_t row = 0; row < kShapeRowCount; ++row) + for (std::uint32_t local = 0; local < kRegularControlCount; + ++local) + { + const std::int32_t sourceId = oneRingSourceIds[ + evaluated * kRegularControlCount + local]; + if (sourceId < 0 || + static_cast(sourceId) >= vertexCount) + { + record_membrane_status(diagnostics, 1, evaluated, sample); + return; + } + const double weight = shapeWeights[ + (sample * kShapeRowCount + row) * kRegularControlCount + + local]; + for (std::uint32_t axis = 0; axis < 3; ++axis) + rows[row][axis] += weight * coordinates[ + static_cast(sourceId) * 3U + axis]; + } + + const double *x = rows[0]; + const double *a_1 = rows[1]; + const double *a_2 = rows[2]; + const double *a_11 = rows[3]; + const double *a_22 = rows[4]; + const double *a_12 = rows[5]; + const double *a_21 = rows[6]; + double xa[3]; + cross3(a_1, a_2, xa); + const double sqa = sqrt(dot3(xa, xa)); + if (!(sqa > 0.0) || !isfinite(sqa)) + { + record_membrane_status(diagnostics, 2, evaluated, sample); + return; + } + const double inverseSqa = 1.0 / sqa; + const double inverseSqaSquared = inverseSqa * inverseSqa; + double tempLeft[3]; + double tempRight[3]; + double xa_1[3]; + double xa_2[3]; + cross3(a_11, a_2, tempLeft); + cross3(a_1, a_21, tempRight); + add3(tempLeft, tempRight, xa_1); + cross3(a_12, a_2, tempLeft); + cross3(a_1, a_22, tempRight); + add3(tempLeft, tempRight, xa_2); + const double sqa_1 = dot3(xa, xa_1) * inverseSqa; + const double sqa_2 = dot3(xa, xa_2) * inverseSqa; + double a_3[3]; + double a_31[3]; + double a_32[3]; + linear3(xa, inverseSqa, xa, 0.0, a_3); + linear3(xa_1, sqa * inverseSqaSquared, + xa, -sqa_1 * inverseSqaSquared, a_31); + linear3(xa_2, sqa * inverseSqaSquared, + xa, -sqa_2 * inverseSqaSquared, a_32); + double a2x3[3]; + double a3x1[3]; + double a1[3]; + double a2[3]; + cross3(a_2, a_3, a2x3); + cross3(a_3, a_1, a3x1); + linear3(a2x3, inverseSqa, a2x3, 0.0, a1); + linear3(a3x1, inverseSqa, a3x1, 0.0, a2); + double a11[3]; + double a12[3]; + double a21[3]; + double a22[3]; + cross3(a_21, a_3, tempLeft); + cross3(a_2, a_31, tempRight); + add3(tempLeft, tempRight, tempLeft); + linear3(tempLeft, sqa * inverseSqaSquared, + a2x3, -sqa_1 * inverseSqaSquared, a11); + cross3(a_22, a_3, tempLeft); + cross3(a_2, a_32, tempRight); + add3(tempLeft, tempRight, tempLeft); + linear3(tempLeft, sqa * inverseSqaSquared, + a2x3, -sqa_2 * inverseSqaSquared, a12); + cross3(a_31, a_1, tempLeft); + cross3(a_3, a_11, tempRight); + add3(tempLeft, tempRight, tempLeft); + linear3(tempLeft, sqa * inverseSqaSquared, + a3x1, -sqa_1 * inverseSqaSquared, a21); + cross3(a_32, a_1, tempLeft); + cross3(a_3, a_12, tempRight); + add3(tempLeft, tempRight, tempLeft); + linear3(tempLeft, sqa * inverseSqaSquared, + a3x1, -sqa_2 * inverseSqaSquared, a22); + const double meanCurvature = + 0.5 * (dot3(a1, a_31) + dot3(a2, a_32)); + const double curvatureDifference = + 2.0 * meanCurvature - spontaneousCurvatures[evaluated]; + const double bendingEnergy = + 0.5 * parameter.kCurv * sqa * curvatureDifference * + curvatureDifference; + const double bendGradientFactor = + -parameter.kCurv * curvatureDifference; + const double bendAreaFactor = + 0.5 * parameter.kCurv * curvatureDifference * + curvatureDifference; + double n1Bend[3]; + double n2Bend[3]; + double m1Bend[3]; + double m2Bend[3]; + for (int axis = 0; axis < 3; ++axis) + { + n1Bend[axis] = bendGradientFactor * + (dot3(a1, a1) * a_31[axis] + + dot3(a1, a2) * a_32[axis]) + + bendAreaFactor * a1[axis]; + n2Bend[axis] = bendGradientFactor * + (dot3(a2, a1) * a_31[axis] + + dot3(a2, a2) * a_32[axis]) + + bendAreaFactor * a2[axis]; + m1Bend[axis] = + parameter.kCurv * curvatureDifference * a1[axis]; + m2Bend[axis] = + parameter.kCurv * curvatureDifference * a2[axis]; + } + double n1Area[3]; + double n2Area[3]; + double n1Volume[3]; + double n2Volume[3]; + for (int axis = 0; axis < 3; ++axis) + { + n1Area[axis] = areaFactor * a1[axis]; + n2Area[axis] = areaFactor * a2[axis]; + n1Volume[axis] = volumeFactor * + (dot3(x, a_3) * a1[axis] - dot3(x, a1) * a_3[axis]); + n2Volume[axis] = volumeFactor * + (dot3(x, a_3) * a2[axis] - dot3(x, a2) * a_3[axis]); + } + if (!finite3(xa_1) || !finite3(xa_2) || !finite3(a_3) || + !finite3(a_31) || !finite3(a_32) || !finite3(a1) || + !finite3(a2) || !finite3(a11) || !finite3(a12) || + !finite3(a21) || !finite3(a22) || + !isfinite(meanCurvature) || !isfinite(bendingEnergy) || + !finite3(n1Bend) || !finite3(n2Bend) || + !finite3(m1Bend) || !finite3(m2Bend) || + !finite3(n1Area) || !finite3(n2Area) || + !finite3(n1Volume) || !finite3(n2Volume)) + { + record_membrane_status(diagnostics, 3, evaluated, sample); + return; + } + + const double coefficient = quadratureCoefficients[sample]; + const double halfCoefficient = 0.5 * coefficient; + area += halfCoefficient * sqa; + volume += kLegacyVolumeQuadratureFactor * coefficient * x[0] * xa[0]; + faceBending += halfCoefficient * bendingEnergy; + faceMean += halfCoefficient * meanCurvature; + for (int axis = 0; axis < 3; ++axis) + accumulatedNormal[axis] += halfCoefficient * a_3[axis]; + const std::uint64_t sampleIndex = + evaluated * kQuadratureSampleCount + sample; + sampleSurfaceMeasures[sampleIndex] = sqa; + sampleMeanCurvatures[sampleIndex] = meanCurvature; + sampleBendingEnergies[sampleIndex] = bendingEnergy; + for (int axis = 0; axis < 3; ++axis) + sampleNormals[sampleIndex * 3U + axis] = a_3[axis]; + + for (std::uint32_t local = 0; local < kRegularControlCount; ++local) + { + const double sf0 = shapeWeights[ + (sample * kShapeRowCount + 0U) * kRegularControlCount + local]; + const double sf1 = shapeWeights[ + (sample * kShapeRowCount + 1U) * kRegularControlCount + local]; + const double sf2 = shapeWeights[ + (sample * kShapeRowCount + 2U) * kRegularControlCount + local]; + const double sf3 = shapeWeights[ + (sample * kShapeRowCount + 3U) * kRegularControlCount + local]; + const double sf4 = shapeWeights[ + (sample * kShapeRowCount + 4U) * kRegularControlCount + local]; + const double sf5 = shapeWeights[ + (sample * kShapeRowCount + 5U) * kRegularControlCount + local]; + const double sf6 = shapeWeights[ + (sample * kShapeRowCount + 6U) * kRegularControlCount + local]; + double da1[3][3]{}; + add_outer_scaled(da1, a1, a_3, -sf3); + add_outer_scaled(da1, a11, a_3, -sf1); + add_outer_scaled(da1, a1, a_31, -sf1); + add_outer_scaled(da1, a2, a_3, -sf6); + add_outer_scaled(da1, a21, a_3, -sf2); + add_outer_scaled(da1, a2, a_31, -sf2); + double da2[3][3]{}; + add_outer_scaled(da2, a1, a_3, -sf5); + add_outer_scaled(da2, a12, a_3, -sf1); + add_outer_scaled(da2, a1, a_32, -sf1); + add_outer_scaled(da2, a2, a_3, -sf4); + add_outer_scaled(da2, a22, a_3, -sf2); + add_outer_scaled(da2, a2, a_32, -sf2); + double da1M1[3]; + double da2M2[3]; + transpose_multiply3(da1, m1Bend, da1M1); + transpose_multiply3(da2, m2Bend, da2M2); + double bending[3]; + double areaForce[3]; + double volumeForce[3]; + for (int axis = 0; axis < 3; ++axis) + { + bending[axis] = -sqa * halfCoefficient * + (da1M1[axis] + da2M2[axis] + + sf1 * n1Bend[axis] + sf2 * n2Bend[axis]); + areaForce[axis] = -sqa * halfCoefficient * + (sf1 * n1Area[axis] + sf2 * n2Area[axis]); + volumeForce[axis] = -sqa * halfCoefficient * + (sf1 * n1Volume[axis] + sf2 * n2Volume[axis] + + volumeFactor * sf0 * a_3[axis]); + } + if (!finite3(bending) || !finite3(areaForce) || + !finite3(volumeForce)) + { + record_membrane_status(diagnostics, 4, evaluated, sample); + return; + } + const std::uint64_t base = + (evaluated * kRegularControlCount + local) * 9U; + for (int axis = 0; axis < 3; ++axis) + { + occurrenceForces[base + axis] += bending[axis]; + occurrenceForces[base + 3U + axis] += areaForce[axis]; + occurrenceForces[base + 6U + axis] += volumeForce[axis]; + } + } + } + const double normalNorm = sqrt(dot3(accumulatedNormal, accumulatedNormal)); + if (!(normalNorm > 0.0) || !isfinite(normalNorm)) + { + record_membrane_status(diagnostics, 2, evaluated, + kQuadratureSampleCount); + return; + } + for (int axis = 0; axis < 3; ++axis) + faceNormals[static_cast(faceId) * 3U + axis] = + accumulatedNormal[axis] / normalNorm; + if (!isfinite(area) || area < 0.0 || !isfinite(volume) || + !isfinite(faceBending) || !isfinite(faceMean)) + { + record_membrane_status(diagnostics, 4, evaluated, + kQuadratureSampleCount); + return; + } + faceAreas[faceId] = area; + faceVolumes[faceId] = volume; + faceBendingEnergies[faceId] = faceBending; + faceMeanCurvatures[faceId] = faceMean; +} + +__global__ void deterministic_membrane_reduction_kernel( + const double *faceAreas, + const double *faceVolumes, + std::int32_t *diagnostics, + double *totals, + std::uint64_t faceCount) +{ + if (blockIdx.x != 0 || threadIdx.x != 0) + return; + double area = 0.0; + double volume = 0.0; + for (std::uint64_t face = 0; face < faceCount; ++face) + { + area += faceAreas[face]; + volume += faceVolumes[face]; + } + if (!isfinite(area) || area < 0.0 || !isfinite(volume)) + record_membrane_status(diagnostics, 4, 0, + kQuadratureSampleCount); + totals[0] = area; + totals[1] = volume; +} + detail::DriverStatus runtime_status(cudaError_t code, const char *operation) { if (code == cudaSuccess) @@ -241,6 +639,97 @@ class RuntimeDriver return runtime_status(cudaGetLastError(), "deterministic_geometry_reduction_kernel"); }; + ops.computeMembrane = [this](const detail::MembraneLaunch &launch) { + cudaError_t code = cudaSetDevice(deviceOrdinal_); + if (code != cudaSuccess) + return runtime_status(code, "cudaSetDevice(compute_membrane)"); + const std::size_t faceBytes = + static_cast(launch.faceCount) * sizeof(double); + const std::size_t evaluated = + static_cast(launch.evaluatedFaceCount); + const std::size_t occurrenceBytes = + evaluated * kRegularControlCount * 9U * sizeof(double); + const std::size_t sampleBytes = + evaluated * kQuadratureSampleCount * sizeof(double); + const std::size_t sampleVectorBytes = sampleBytes * 3U; + const auto clear = [this](detail::DeviceBufferHandle handle, + std::size_t bytes) { + if (bytes == 0) + return cudaSuccess; + return cudaMemsetAsync(reinterpret_cast(handle), 0, + bytes, stream_); + }; + code = clear(launch.faceAreas, faceBytes); + if (code == cudaSuccess) + code = clear(launch.faceVolumes, faceBytes); + if (code == cudaSuccess) + code = clear(launch.geometryTotals, 2U * sizeof(double)); + if (code == cudaSuccess) + code = clear(launch.faceBendingEnergies, faceBytes); + if (code == cudaSuccess) + code = clear(launch.faceMeanCurvatures, faceBytes); + if (code == cudaSuccess) + code = clear(launch.faceNormals, faceBytes * 3U); + if (code == cudaSuccess) + code = clear(launch.occurrenceForces, occurrenceBytes); + if (code == cudaSuccess) + code = clear(launch.sampleSurfaceMeasures, sampleBytes); + if (code == cudaSuccess) + code = clear(launch.sampleMeanCurvatures, sampleBytes); + if (code == cudaSuccess) + code = clear(launch.sampleNormals, sampleVectorBytes); + if (code == cudaSuccess) + code = clear(launch.sampleBendingEnergies, sampleBytes); + if (code == cudaSuccess) + code = clear(launch.statusDiagnostics, + 3U * sizeof(std::int32_t)); + if (code != cudaSuccess) + return runtime_status(code, "cudaMemsetAsync(membrane)"); + constexpr unsigned int threads = 128; + const unsigned int blocks = static_cast( + (launch.evaluatedFaceCount + threads - 1) / threads); + if (blocks != 0) + { + regular_membrane_kernel<<>>( + reinterpret_cast( + launch.evaluatedFaceIds), + reinterpret_cast( + launch.oneRingSourceIds), + reinterpret_cast( + launch.spontaneousCurvatures), + reinterpret_cast( + launch.packedParameters), + reinterpret_cast( + launch.quadratureCoefficients), + reinterpret_cast(launch.shapeWeights), + reinterpret_cast(launch.coordinates), + reinterpret_cast(launch.faceAreas), + reinterpret_cast(launch.faceVolumes), + reinterpret_cast(launch.faceBendingEnergies), + reinterpret_cast(launch.faceMeanCurvatures), + reinterpret_cast(launch.faceNormals), + reinterpret_cast(launch.occurrenceForces), + reinterpret_cast(launch.sampleSurfaceMeasures), + reinterpret_cast(launch.sampleMeanCurvatures), + reinterpret_cast(launch.sampleNormals), + reinterpret_cast(launch.sampleBendingEnergies), + reinterpret_cast( + launch.statusDiagnostics), + launch.vertexCount, launch.faceCount, + launch.evaluatedFaceCount); + code = cudaGetLastError(); + if (code != cudaSuccess) + return runtime_status(code, "regular_membrane_kernel"); + } + deterministic_membrane_reduction_kernel<<<1, 1, 0, stream_>>>( + reinterpret_cast(launch.faceAreas), + reinterpret_cast(launch.faceVolumes), + reinterpret_cast(launch.statusDiagnostics), + reinterpret_cast(launch.geometryTotals), + launch.faceCount); + return runtime_status(cudaGetLastError(), + "deterministic_membrane_reduction_kernel"); + }; ops.synchronize = [this]() { return runtime_status(cudaStreamSynchronize(stream_), "cudaStreamSynchronize"); diff --git a/src/cuda/Cuda_mesh_state_common.cpp b/src/cuda/Cuda_mesh_state_common.cpp index 3f48d7a..de28965 100644 --- a/src/cuda/Cuda_mesh_state_common.cpp +++ b/src/cuda/Cuda_mesh_state_common.cpp @@ -46,6 +46,8 @@ const char *transfer_reason_name(TransferReason reason) noexcept case TransferReason::CandidateCoordinates: return "candidate_coordinates"; case TransferReason::CandidateGeometry: return "candidate_geometry"; case TransferReason::GeometryDiagnostics: return "geometry_diagnostics"; + case TransferReason::CandidateMembrane: return "candidate_membrane"; + case TransferReason::MembraneDiagnostics: return "membrane_diagnostics"; case TransferReason::Count: break; } return "unknown"; @@ -280,6 +282,7 @@ struct MeshStateCore::Impl BufferGroup coordinates; BufferGroup reference; BufferGroup geometry; + BufferGroup membrane; BufferGroup deferredCleanup; std::uint64_t vertexCount = 0; std::uint64_t faceCount = 0; @@ -401,6 +404,49 @@ struct MeshStateCore::Impl return result; } + std::vector membrane_views(const RegularMeshPack &pack, + bool &ok) const + { + std::vector result(9); + std::uint64_t faceDoubles = 0; + std::uint64_t occurrenceCount = 0; + std::uint64_t occurrenceDoubles = 0; + std::uint64_t sampleCount = 0; + std::uint64_t sampleVectorDoubles = 0; + ok = multiply(pack.faceCount, sizeof(double), faceDoubles) && + multiply(pack.evaluatedFaceCount, kRegularControlCount, + occurrenceCount) && + multiply(occurrenceCount, 9U, occurrenceDoubles) && + multiply(occurrenceDoubles, sizeof(double), occurrenceDoubles) && + multiply(pack.evaluatedFaceCount, kQuadratureSampleCount, + sampleCount) && + multiply(sampleCount, 3U, sampleVectorDoubles) && + multiply(sampleVectorDoubles, sizeof(double), + sampleVectorDoubles) && + faceDoubles <= std::numeric_limits::max() && + occurrenceDoubles <= std::numeric_limits::max() && + sampleCount <= std::numeric_limits::max() / + sizeof(double) && + sampleVectorDoubles <= std::numeric_limits::max(); + ok = ok && faceDoubles <= + std::numeric_limits::max() / 3U; + if (!ok) + return result; + const std::size_t faces = static_cast(faceDoubles); + const std::size_t samples = + static_cast(sampleCount) * sizeof(double); + result[0].bytes = faces; + result[1].bytes = faces; + result[2].bytes = faces * 3U; + result[3].bytes = static_cast(occurrenceDoubles); + result[4].bytes = samples; + result[5].bytes = samples; + result[6].bytes = static_cast(sampleVectorDoubles); + result[7].bytes = samples; + result[8].bytes = 3U * sizeof(std::int32_t); + return result; + } + DeviceStateError stage_group(PendingGroup &pending) { const BufferGroup &old = *pending.destination; @@ -583,7 +629,8 @@ struct MeshStateCore::Impl { report.residentBytes = group_bytes(topology) + group_bytes(numericalPlan) + group_bytes(parameters) + group_bytes(coordinates) + - group_bytes(reference) + group_bytes(geometry); + group_bytes(reference) + group_bytes(geometry) + + group_bytes(membrane); report.capacityBytes[static_cast(TransferReason::Topology)] = group_bytes(topology); report.capacityBytes[static_cast( @@ -601,6 +648,8 @@ struct MeshStateCore::Impl : 0; report.capacityBytes[static_cast( TransferReason::CandidateGeometry)] = group_bytes(geometry); + report.capacityBytes[static_cast( + TransferReason::CandidateMembrane)] = group_bytes(membrane); } }; @@ -683,6 +732,10 @@ DeviceStateError MeshStateCore::ensure_resident(const RegularMeshPack &pack) pending.push_back({&s.geometry, {}, TransferReason::CandidateGeometry, s.geometry_views(pack, viewsOk)}); + if (topologyChanged) + pending.push_back({&s.membrane, {}, + TransferReason::CandidateMembrane, + s.membrane_views(pack, viewsOk)}); if (!viewsOk) return s.record(error(DeviceStateErrorCode::ArithmeticOverflow, "ensure_resident", "host buffer byte count overflowed")); @@ -931,6 +984,266 @@ GeometryCandidateResult MeshStateCore::compute_candidate_geometry() return result; } +MembraneCandidateResult MeshStateCore::compute_candidate_membrane() +{ + auto &s = *impl_; + MembraneCandidateResult result; + result.coordinateGeneration = s.report.candidateGeneration; + s.report.lastDirtyGroups.fill(false); + if (s.report.phase != TransactionPhase::CandidatePrepared || + s.topology.buffers.size() < 7 || s.parameters.buffers.size() < 3 || + s.numericalPlan.buffers.size() < 3 || + s.coordinates.buffers.size() != 3 || s.geometry.buffers.size() != 4 || + s.membrane.buffers.size() != 9) + { + result.error = error(DeviceStateErrorCode::InvalidTransition, + "compute_candidate_membrane", + "membrane evaluation requires a prepared candidate and complete resident buffers"); + s.record(result.error); + return result; + } + if (!s.operations.computeMembrane) + { + result.error = error(DeviceStateErrorCode::InvalidConfiguration, + "compute_candidate_membrane", + "the device driver does not provide the Step-5 membrane operation"); + s.record(result.error); + return result; + } + if (s.report.cleanupPending) + { + result.error = s.cleanup_blocked("compute_candidate_membrane"); + s.record(result.error); + return result; + } + + s.report.phase = TransactionPhase::Computing; + s.report.lastDirtyGroups[static_cast( + TransferReason::CandidateMembrane)] = true; + const MembraneLaunch launch{ + s.topology.buffers[4].handle, + s.topology.buffers[6].handle, + s.parameters.buffers[1].handle, + s.parameters.buffers[2].handle, + s.numericalPlan.buffers[1].handle, + s.numericalPlan.buffers[2].handle, + s.coordinates.buffers[s.report.candidateCoordinateSlot].handle, + s.geometry.buffers[0].handle, + s.geometry.buffers[1].handle, + s.geometry.buffers[3].handle, + s.membrane.buffers[0].handle, + s.membrane.buffers[1].handle, + s.membrane.buffers[2].handle, + s.membrane.buffers[3].handle, + s.membrane.buffers[4].handle, + s.membrane.buffers[5].handle, + s.membrane.buffers[6].handle, + s.membrane.buffers[7].handle, + s.membrane.buffers[8].handle, + s.vertexCount, + s.faceCount, + s.evaluatedFaceCount, + }; + DriverStatus driverStatus = s.operations.computeMembrane(launch); + if (!driverStatus.success) + { + result.error = error(DeviceStateErrorCode::CandidateFailed, + driverStatus.operation.empty() + ? "compute_membrane" + : driverStatus.operation, + driverStatus.message, driverStatus.nativeCode); + s.report.phase = TransactionPhase::Failed; + s.report.lastOutcome = TransactionOutcome::Failed; + s.record(result.error); + return result; + } + driverStatus = s.operations.synchronize(); + if (!driverStatus.success) + { + result.error = error(DeviceStateErrorCode::SynchronizationFailed, + driverStatus.operation.empty() + ? "synchronize_membrane" + : driverStatus.operation, + driverStatus.message, driverStatus.nativeCode); + s.report.phase = TransactionPhase::Failed; + s.report.lastOutcome = TransactionOutcome::Failed; + s.record(result.error); + return result; + } + ++s.report.synchronizations; + + const std::size_t faces = static_cast(s.faceCount); + const std::size_t evaluated = + static_cast(s.evaluatedFaceCount); + result.faceAreas.resize(faces); + result.faceVolumes.resize(faces); + result.faceBendingEnergies.resize(faces); + result.faceMeanCurvatures.resize(faces); + result.faceNormals.resize(faces * 3U); + result.occurrenceForces.resize(evaluated * kRegularControlCount * 9U); + result.sampleSurfaceMeasures.resize(evaluated * kQuadratureSampleCount); + result.sampleMeanCurvatures.resize(evaluated * kQuadratureSampleCount); + result.sampleNormals.resize(evaluated * kQuadratureSampleCount * 3U); + result.sampleBendingEnergies.resize(evaluated * kQuadratureSampleCount); + double totals[2]{}; + std::int32_t diagnostics[3]{}; + auto &counter = s.report.transfers[static_cast( + TransferReason::MembraneDiagnostics)]; + std::uint64_t completedOperations = 0; + std::uint64_t diagnosticBytes = 0; + const auto copy = [&](void *destination, DeviceBufferHandle source, + std::size_t bytes) -> bool { + if (bytes == 0) + return true; + ++counter.attemptedOperations; + counter.attemptedBytes += bytes; + const DriverStatus copied = + s.operations.copyDeviceToHost(destination, source, bytes); + if (!copied.success) + { + result.error = error( + DeviceStateErrorCode::TransferFailed, + copied.operation.empty() ? "copy_membrane_diagnostic" + : copied.operation, + copied.message, copied.nativeCode); + return false; + } + ++completedOperations; + diagnosticBytes += bytes; + return true; + }; + const std::size_t faceBytes = faces * sizeof(double); + if (!copy(result.faceAreas.data(), s.geometry.buffers[0].handle, + faceBytes) || + !copy(result.faceVolumes.data(), s.geometry.buffers[1].handle, + faceBytes) || + !copy(totals, s.geometry.buffers[3].handle, sizeof(totals)) || + !copy(result.faceBendingEnergies.data(), s.membrane.buffers[0].handle, + faceBytes) || + !copy(result.faceMeanCurvatures.data(), s.membrane.buffers[1].handle, + faceBytes) || + !copy(result.faceNormals.data(), s.membrane.buffers[2].handle, + result.faceNormals.size() * sizeof(double)) || + !copy(result.occurrenceForces.data(), s.membrane.buffers[3].handle, + result.occurrenceForces.size() * sizeof(double)) || + !copy(result.sampleSurfaceMeasures.data(), + s.membrane.buffers[4].handle, + result.sampleSurfaceMeasures.size() * sizeof(double)) || + !copy(result.sampleMeanCurvatures.data(), + s.membrane.buffers[5].handle, + result.sampleMeanCurvatures.size() * sizeof(double)) || + !copy(result.sampleNormals.data(), s.membrane.buffers[6].handle, + result.sampleNormals.size() * sizeof(double)) || + !copy(result.sampleBendingEnergies.data(), + s.membrane.buffers[7].handle, + result.sampleBendingEnergies.size() * sizeof(double)) || + !copy(diagnostics, s.membrane.buffers[8].handle, + sizeof(diagnostics))) + { + s.report.phase = TransactionPhase::Failed; + s.report.lastOutcome = TransactionOutcome::Failed; + s.record(result.error); + return result; + } + driverStatus = s.operations.synchronize(); + if (!driverStatus.success) + { + result.error = error(DeviceStateErrorCode::SynchronizationFailed, + driverStatus.operation.empty() + ? "synchronize_membrane_diagnostic" + : driverStatus.operation, + driverStatus.message, driverStatus.nativeCode); + s.report.phase = TransactionPhase::Failed; + s.report.lastOutcome = TransactionOutcome::Failed; + s.record(result.error); + return result; + } + ++s.report.synchronizations; + counter.completedOperations += completedOperations; + counter.completedBytes += diagnosticBytes; + result.totalArea = totals[0]; + result.totalVolume = totals[1]; + result.status = static_cast(diagnostics[0]); + result.failedEvaluatedFace = diagnostics[1] < 0 + ? 0U + : static_cast(diagnostics[1]); + result.failedSample = diagnostics[2] < 0 + ? 0U + : static_cast(diagnostics[2]); + const auto allFinite = [](const std::vector &values) { + return std::all_of(values.begin(), values.end(), + [](double value) { return std::isfinite(value); }); + }; + const auto normalsValid = [](const std::vector &values, + bool allowZero) { + if (values.size() % 3U != 0U) + return false; + for (std::size_t index = 0; index < values.size(); index += 3U) + { + const double squared = values[index] * values[index] + + values[index + 1U] * values[index + 1U] + + values[index + 2U] * values[index + 2U]; + if (allowZero && squared == 0.0) + continue; + if (!std::isfinite(squared) || std::abs(squared - 1.0) > 1.0e-10) + return false; + } + return true; + }; + double summedArea = 0.0; + double summedVolume = 0.0; + for (std::size_t face = 0; face < result.faceAreas.size(); ++face) + { + summedArea += result.faceAreas[face]; + summedVolume += result.faceVolumes[face]; + } + const auto totalMatches = [](double actual, double expected) { + return std::abs(actual - expected) <= + 1.0e-12 * std::max(1.0, std::abs(expected)); + }; + const bool valid = + result.status == MembraneStatusCode::None && + allFinite(result.faceAreas) && + std::all_of(result.faceAreas.begin(), result.faceAreas.end(), + [](double value) { return value >= 0.0; }) && + allFinite(result.faceVolumes) && + allFinite(result.faceBendingEnergies) && + std::all_of(result.faceBendingEnergies.begin(), + result.faceBendingEnergies.end(), + [](double value) { return value >= 0.0; }) && + allFinite(result.faceMeanCurvatures) && allFinite(result.faceNormals) && + normalsValid(result.faceNormals, true) && + allFinite(result.occurrenceForces) && + allFinite(result.sampleSurfaceMeasures) && + std::all_of(result.sampleSurfaceMeasures.begin(), + result.sampleSurfaceMeasures.end(), + [](double value) { return value > 0.0; }) && + allFinite(result.sampleMeanCurvatures) && + allFinite(result.sampleNormals) && + normalsValid(result.sampleNormals, false) && + allFinite(result.sampleBendingEnergies) && + std::all_of(result.sampleBendingEnergies.begin(), + result.sampleBendingEnergies.end(), + [](double value) { return value >= 0.0; }) && + std::isfinite(result.totalArea) && result.totalArea >= 0.0 && + std::isfinite(result.totalVolume) && + totalMatches(result.totalArea, summedArea) && + totalMatches(result.totalVolume, summedVolume); + if (!valid) + { + result.error = error(DeviceStateErrorCode::CandidateFailed, + "validate_candidate_membrane", + "device membrane evaluation reported a degeneracy or nonfinite output"); + s.report.phase = TransactionPhase::Failed; + s.report.lastOutcome = TransactionOutcome::Failed; + s.record(result.error); + return result; + } + s.report.phase = TransactionPhase::Validated; + s.clear_error(); + return result; +} + DeviceStateError MeshStateCore::mark_computing() { auto &s = *impl_; @@ -1076,7 +1389,7 @@ DeviceStateError MeshStateCore::close() auto &s = *impl_; if (!s.closing) { - for (BufferGroup *group : {&s.geometry, &s.reference, &s.coordinates, &s.parameters, + for (BufferGroup *group : {&s.membrane, &s.geometry, &s.reference, &s.coordinates, &s.parameters, &s.numericalPlan, &s.topology}) { for (auto &buffer : group->buffers) @@ -1241,6 +1554,20 @@ GeometryCandidateResult CudaMeshState::compute_candidate_geometry() return result; } +MembraneCandidateResult CudaMeshState::compute_candidate_membrane() +{ + MembraneCandidateResult result; + DeviceStateError blocked = impl_->guard("compute_candidate_membrane"); + if (!blocked.ok()) + { + result.error = blocked; + return result; + } + result = impl_->core->compute_candidate_membrane(); + impl_->refresh(); + return result; +} + DeviceStateError CudaMeshState::mark_computing() { DeviceStateError blocked = impl_->guard("mark_computing"); diff --git a/src/cuda/Cuda_regular_membrane_cpu.cpp b/src/cuda/Cuda_regular_membrane_cpu.cpp new file mode 100644 index 0000000..97e4bd2 --- /dev/null +++ b/src/cuda/Cuda_regular_membrane_cpu.cpp @@ -0,0 +1,397 @@ +#include "cuda/detail/Cuda_regular_membrane_cpu.hpp" + +#include +#include +#include +#include +#include + +namespace slimed::cuda_residency::detail +{ +namespace +{ +using Vec3 = std::array; +using Mat3 = std::array; + +constexpr double kLegacyVolumeQuadratureFactor = 0.16666666666; + +Vec3 add(const Vec3 &left, const Vec3 &right) +{ + return {left[0] + right[0], left[1] + right[1], left[2] + right[2]}; +} + +Vec3 subtract(const Vec3 &left, const Vec3 &right) +{ + return {left[0] - right[0], left[1] - right[1], left[2] - right[2]}; +} + +Vec3 scale(const Vec3 &value, double factor) +{ + return {value[0] * factor, value[1] * factor, value[2] * factor}; +} + +double dot(const Vec3 &left, const Vec3 &right) +{ + return left[0] * right[0] + left[1] * right[1] + + left[2] * right[2]; +} + +Vec3 cross(const Vec3 &left, const Vec3 &right) +{ + return {left[1] * right[2] - left[2] * right[1], + left[2] * right[0] - left[0] * right[2], + left[0] * right[1] - left[1] * right[0]}; +} + +Mat3 outer(const Vec3 &left, const Vec3 &right) +{ + Mat3 result{}; + for (std::size_t row = 0; row < 3; ++row) + for (std::size_t column = 0; column < 3; ++column) + result[row][column] = left[row] * right[column]; + return result; +} + +void add_scaled(Mat3 &destination, const Mat3 &source, double factor) +{ + for (std::size_t row = 0; row < 3; ++row) + for (std::size_t column = 0; column < 3; ++column) + destination[row][column] += factor * source[row][column]; +} + +Vec3 transpose_multiply(const Mat3 &matrix, const Vec3 &value) +{ + Vec3 result{}; + for (std::size_t column = 0; column < 3; ++column) + for (std::size_t row = 0; row < 3; ++row) + result[column] += value[row] * matrix[row][column]; + return result; +} + +bool finite(const Vec3 &value) +{ + return std::isfinite(value[0]) && std::isfinite(value[1]) && + std::isfinite(value[2]); +} + +bool finite(const Mat3 &value) +{ + return finite(value[0]) && finite(value[1]) && finite(value[2]); +} + +RegularMembraneCpuResult fail(RegularMembraneCpuResult result, + RegularMembraneStatus status, + std::size_t evaluated, + std::size_t sample, + const char *message) +{ + result.status = status; + result.failedEvaluatedFace = evaluated; + result.failedSample = static_cast(sample); + result.message = message; + return result; +} + +} // namespace + +const char *regular_membrane_status_name(RegularMembraneStatus status) noexcept +{ + switch (status) + { + case RegularMembraneStatus::None: return "none"; + case RegularMembraneStatus::InvalidSource: return "invalid_source"; + case RegularMembraneStatus::DegenerateSample: return "degenerate_sample"; + case RegularMembraneStatus::NonFiniteIntermediate: + return "nonfinite_intermediate"; + case RegularMembraneStatus::NonFiniteOutput: return "nonfinite_output"; + } + return "unknown"; +} + +RegularMembraneCpuResult evaluate_regular_membrane_cpu( + const RegularMeshPack &pack, + const std::vector &coordinates) +{ + RegularMembraneCpuResult result; + const std::size_t faceCount = static_cast(pack.faceCount); + const std::size_t evaluatedCount = + static_cast(pack.evaluatedFaceCount); + result.faceAreas.assign(faceCount, 0.0); + result.faceVolumes.assign(faceCount, 0.0); + result.faceBendingEnergies.assign(faceCount, 0.0); + result.faceMeanCurvatures.assign(faceCount, 0.0); + result.faceNormals.assign(faceCount * 3, 0.0); + result.occurrenceForces.assign( + evaluatedCount * kRegularControlCount * 9U, 0.0); + result.sampleSurfaceMeasures.assign( + evaluatedCount * kQuadratureSampleCount, 0.0); + result.sampleMeanCurvatures.assign( + evaluatedCount * kQuadratureSampleCount, 0.0); + result.sampleNormals.assign( + evaluatedCount * kQuadratureSampleCount * 3U, 0.0); + result.sampleBendingEnergies.assign( + evaluatedCount * kQuadratureSampleCount, 0.0); + + const double uSurfPerArea = + pack.parameters.uSurf == 0.0 || pack.parameters.area0 == 0.0 + ? 0.0 + : pack.parameters.uSurf / pack.parameters.area0; + const double uVol = + pack.parameters.uVol == 0.0 || pack.parameters.vol0 == 0.0 + ? 0.0 + : pack.parameters.uVol / pack.parameters.vol0; + const double areaFactor = + uSurfPerArea * (pack.parameters.area - pack.parameters.area0); + const double volumeFactor = + uVol * (pack.parameters.vol - pack.parameters.vol0) / 3.0; + + for (std::size_t evaluated = 0; evaluated < evaluatedCount; ++evaluated) + { + const std::int32_t faceId = pack.evaluatedFaceIds[evaluated]; + if (faceId < 0 || static_cast(faceId) >= faceCount) + return fail(std::move(result), + RegularMembraneStatus::InvalidSource, evaluated, 0, + "evaluated face ID is outside the declared face range"); + const std::size_t face = static_cast(faceId); + Vec3 accumulatedNormal{}; + for (std::size_t sample = 0; sample < kQuadratureSampleCount; ++sample) + { + Vec3 rows[kShapeRowCount]{}; + for (std::size_t row = 0; row < kShapeRowCount; ++row) + for (std::size_t local = 0; local < kRegularControlCount; + ++local) + { + const std::int32_t sourceId = pack.oneRingSourceIds[ + evaluated * kRegularControlCount + local]; + if (sourceId < 0 || + static_cast(sourceId) >= + static_cast(pack.vertexCount)) + return fail(std::move(result), + RegularMembraneStatus::InvalidSource, + evaluated, sample, + "one-ring source ID is outside the vertex range"); + const double weight = pack.shapeWeights[ + (sample * kShapeRowCount + row) * + kRegularControlCount + + local]; + const std::size_t source = + static_cast(sourceId); + for (std::size_t axis = 0; axis < 3; ++axis) + rows[row][axis] += + weight * coordinates[source * 3 + axis]; + } + + const Vec3 &x = rows[0]; + const Vec3 &a_1 = rows[1]; + const Vec3 &a_2 = rows[2]; + const Vec3 &a_11 = rows[3]; + const Vec3 &a_22 = rows[4]; + const Vec3 &a_12 = rows[5]; + const Vec3 &a_21 = rows[6]; + const Vec3 xa = cross(a_1, a_2); + const double sqa = std::sqrt(dot(xa, xa)); + if (!(sqa > 0.0) || !std::isfinite(sqa)) + return fail(std::move(result), + RegularMembraneStatus::DegenerateSample, + evaluated, sample, + "regular membrane sample has zero or nonfinite surface measure"); + const double inverseSqa = 1.0 / sqa; + const double inverseSqaSquared = inverseSqa * inverseSqa; + const Vec3 xa_1 = add(cross(a_11, a_2), cross(a_1, a_21)); + const Vec3 xa_2 = add(cross(a_12, a_2), cross(a_1, a_22)); + const double sqa_1 = dot(xa, xa_1) * inverseSqa; + const double sqa_2 = dot(xa, xa_2) * inverseSqa; + const Vec3 a_3 = scale(xa, inverseSqa); + const Vec3 a_31 = scale( + subtract(scale(xa_1, sqa), scale(xa, sqa_1)), + inverseSqaSquared); + const Vec3 a_32 = scale( + subtract(scale(xa_2, sqa), scale(xa, sqa_2)), + inverseSqaSquared); + const Vec3 a2x3 = cross(a_2, a_3); + const Vec3 a3x1 = cross(a_3, a_1); + const Vec3 a1 = scale(a2x3, inverseSqa); + const Vec3 a2 = scale(a3x1, inverseSqa); + const Vec3 a11 = scale( + subtract(scale(add(cross(a_21, a_3), cross(a_2, a_31)), + sqa), + scale(a2x3, sqa_1)), + inverseSqaSquared); + const Vec3 a12 = scale( + subtract(scale(add(cross(a_22, a_3), cross(a_2, a_32)), + sqa), + scale(a2x3, sqa_2)), + inverseSqaSquared); + const Vec3 a21 = scale( + subtract(scale(add(cross(a_31, a_1), cross(a_3, a_11)), + sqa), + scale(a3x1, sqa_1)), + inverseSqaSquared); + const Vec3 a22 = scale( + subtract(scale(add(cross(a_32, a_1), cross(a_3, a_12)), + sqa), + scale(a3x1, sqa_2)), + inverseSqaSquared); + const double meanCurvature = + 0.5 * (dot(a1, a_31) + dot(a2, a_32)); + const double curvatureDifference = + 2.0 * meanCurvature - + pack.evaluatedFaceSpontaneousCurvature[evaluated]; + const double bendingEnergy = + 0.5 * pack.parameters.kCurv * sqa * + curvatureDifference * curvatureDifference; + + const double bendGradientFactor = + -pack.parameters.kCurv * curvatureDifference; + const double bendAreaFactor = + 0.5 * pack.parameters.kCurv * curvatureDifference * + curvatureDifference; + const Vec3 n1Bend = add( + scale(add(scale(a_31, dot(a1, a1)), + scale(a_32, dot(a1, a2))), + bendGradientFactor), + scale(a1, bendAreaFactor)); + const Vec3 n2Bend = add( + scale(add(scale(a_31, dot(a2, a1)), + scale(a_32, dot(a2, a2))), + bendGradientFactor), + scale(a2, bendAreaFactor)); + const Vec3 m1Bend = + scale(a1, pack.parameters.kCurv * curvatureDifference); + const Vec3 m2Bend = + scale(a2, pack.parameters.kCurv * curvatureDifference); + const Vec3 n1Area = scale(a1, areaFactor); + const Vec3 n2Area = scale(a2, areaFactor); + const Vec3 n1Volume = scale( + subtract(scale(a1, dot(x, a_3)), + scale(a_3, dot(x, a1))), + volumeFactor); + const Vec3 n2Volume = scale( + subtract(scale(a2, dot(x, a_3)), + scale(a_3, dot(x, a2))), + volumeFactor); + + if (!finite(xa_1) || !finite(xa_2) || !finite(a_3) || + !finite(a_31) || !finite(a_32) || !finite(a1) || + !finite(a2) || !finite(a11) || !finite(a12) || + !finite(a21) || !finite(a22) || + !std::isfinite(meanCurvature) || + !std::isfinite(bendingEnergy) || !finite(n1Bend) || + !finite(n2Bend) || !finite(m1Bend) || !finite(m2Bend) || + !finite(n1Area) || !finite(n2Area) || !finite(n1Volume) || + !finite(n2Volume)) + return fail(std::move(result), + RegularMembraneStatus::NonFiniteIntermediate, + evaluated, sample, + "regular membrane intermediate is nonfinite"); + + const double coefficient = pack.quadratureCoefficients[sample]; + const double halfCoefficient = 0.5 * coefficient; + result.faceAreas[face] += halfCoefficient * sqa; + result.faceVolumes[face] += + kLegacyVolumeQuadratureFactor * coefficient * x[0] * xa[0]; + result.faceBendingEnergies[face] += + halfCoefficient * bendingEnergy; + result.faceMeanCurvatures[face] += + halfCoefficient * meanCurvature; + accumulatedNormal = + add(accumulatedNormal, scale(a_3, halfCoefficient)); + const std::size_t sampleIndex = + evaluated * kQuadratureSampleCount + sample; + result.sampleSurfaceMeasures[sampleIndex] = sqa; + result.sampleMeanCurvatures[sampleIndex] = meanCurvature; + result.sampleBendingEnergies[sampleIndex] = bendingEnergy; + for (std::size_t axis = 0; axis < 3; ++axis) + result.sampleNormals[sampleIndex * 3 + axis] = a_3[axis]; + + for (std::size_t local = 0; local < kRegularControlCount; + ++local) + { + const double *weights = &pack.shapeWeights[ + sample * kShapeRowCount * kRegularControlCount + local]; + const double sf0 = weights[0 * kRegularControlCount]; + const double sf1 = weights[1 * kRegularControlCount]; + const double sf2 = weights[2 * kRegularControlCount]; + const double sf3 = weights[3 * kRegularControlCount]; + const double sf4 = weights[4 * kRegularControlCount]; + const double sf5 = weights[5 * kRegularControlCount]; + const double sf6 = weights[6 * kRegularControlCount]; + Mat3 da1{}; + add_scaled(da1, outer(a1, a_3), -sf3); + add_scaled(da1, outer(a11, a_3), -sf1); + add_scaled(da1, outer(a1, a_31), -sf1); + add_scaled(da1, outer(a2, a_3), -sf6); + add_scaled(da1, outer(a21, a_3), -sf2); + add_scaled(da1, outer(a2, a_31), -sf2); + Mat3 da2{}; + add_scaled(da2, outer(a1, a_3), -sf5); + add_scaled(da2, outer(a12, a_3), -sf1); + add_scaled(da2, outer(a1, a_32), -sf1); + add_scaled(da2, outer(a2, a_3), -sf4); + add_scaled(da2, outer(a22, a_3), -sf2); + add_scaled(da2, outer(a2, a_32), -sf2); + Vec3 bending = add( + add(transpose_multiply(da1, m1Bend), + transpose_multiply(da2, m2Bend)), + add(scale(n1Bend, sf1), scale(n2Bend, sf2))); + bending = scale(bending, -sqa * halfCoefficient); + Vec3 area = scale( + add(scale(n1Area, sf1), scale(n2Area, sf2)), + -sqa * halfCoefficient); + Vec3 volume = add( + add(scale(n1Volume, sf1), scale(n2Volume, sf2)), + scale(a_3, volumeFactor * sf0)); + volume = scale(volume, -sqa * halfCoefficient); + if (!finite(da1) || !finite(da2) || !finite(bending) || + !finite(area) || !finite(volume)) + return fail(std::move(result), + RegularMembraneStatus::NonFiniteOutput, + evaluated, sample, + "regular membrane occurrence force is nonfinite"); + const std::size_t base = + (evaluated * kRegularControlCount + local) * 9U; + for (std::size_t axis = 0; axis < 3; ++axis) + { + result.occurrenceForces[base + axis] += bending[axis]; + result.occurrenceForces[base + 3 + axis] += area[axis]; + result.occurrenceForces[base + 6 + axis] += volume[axis]; + } + } + } + const double normalNorm = std::sqrt(dot(accumulatedNormal, + accumulatedNormal)); + if (!(normalNorm > 0.0) || !std::isfinite(normalNorm)) + return fail(std::move(result), + RegularMembraneStatus::DegenerateSample, evaluated, + kQuadratureSampleCount, + "integrated face normal is zero or nonfinite"); + accumulatedNormal = scale(accumulatedNormal, 1.0 / normalNorm); + for (std::size_t axis = 0; axis < 3; ++axis) + result.faceNormals[face * 3 + axis] = accumulatedNormal[axis]; + } + for (std::size_t face = 0; face < faceCount; ++face) + { + result.totalArea += result.faceAreas[face]; + result.totalVolume += result.faceVolumes[face]; + } + const auto allFinite = [](const std::vector &values) { + return std::all_of(values.begin(), values.end(), + [](double value) { return std::isfinite(value); }); + }; + if (!allFinite(result.faceAreas) || !allFinite(result.faceVolumes) || + !allFinite(result.faceBendingEnergies) || + !allFinite(result.faceMeanCurvatures) || + !allFinite(result.faceNormals) || !allFinite(result.occurrenceForces) || + !allFinite(result.sampleSurfaceMeasures) || + !allFinite(result.sampleMeanCurvatures) || + !allFinite(result.sampleNormals) || + !allFinite(result.sampleBendingEnergies) || + !std::isfinite(result.totalArea) || + !std::isfinite(result.totalVolume)) + return fail(std::move(result), + RegularMembraneStatus::NonFiniteOutput, 0, 0, + "regular membrane result contains a nonfinite value"); + return result; +} + +} // namespace slimed::cuda_residency::detail diff --git a/tests/test_cuda_mesh_state.cpp b/tests/test_cuda_mesh_state.cpp index 4c0608b..e50675c 100644 --- a/tests/test_cuda_mesh_state.cpp +++ b/tests/test_cuda_mesh_state.cpp @@ -102,6 +102,16 @@ struct FakeDevice NonFiniteTotals, }; + enum class MembraneCorruption + { + None, + Status, + NonFiniteOccurrence, + ZeroSurfaceMeasure, + NonUnitNormal, + MismatchedTotals, + }; + static constexpr int kInjected = 73; std::unordered_map> memory; DeviceBufferHandle nextHandle = 1; @@ -111,15 +121,18 @@ struct FakeDevice std::uint64_t copyCalls = 0; std::uint64_t copyToHostCalls = 0; std::uint64_t geometryCalls = 0; + std::uint64_t membraneCalls = 0; std::uint64_t synchronizeCalls = 0; std::uint64_t releaseCalls = 0; std::uint64_t failAllocationCall = 0; std::uint64_t failCopyCall = 0; std::uint64_t failCopyToHostCall = 0; std::uint64_t failGeometryCall = 0; + std::uint64_t failMembraneCall = 0; std::uint64_t failSynchronizeCall = 0; std::uint64_t failReleaseCall = 0; GeometryCorruption geometryCorruption = GeometryCorruption::None; + MembraneCorruption membraneCorruption = MembraneCorruption::None; DeviceOperations operations() { @@ -283,6 +296,80 @@ struct FakeDevice } return DriverStatus{}; }; + ops.computeMembrane = [this](const MembraneLaunch &launch) { + ++membraneCalls; + if (membraneCalls == failMembraneCall) + return DriverStatus{false, kInjected, + "fake_compute_membrane", + "injected membrane failure"}; + const auto clear = [this](DeviceBufferHandle handle) { + std::fill(memory.at(handle).begin(), memory.at(handle).end(), + 0); + }; + const auto write_double = [this](DeviceBufferHandle handle, + std::size_t index, + double value) { + std::memcpy(memory.at(handle).data() + index * sizeof(value), + &value, sizeof(value)); + }; + const auto write_i32 = [this](DeviceBufferHandle handle, + std::size_t index, + std::int32_t value) { + std::memcpy(memory.at(handle).data() + index * sizeof(value), + &value, sizeof(value)); + }; + for (const DeviceBufferHandle handle : { + launch.faceAreas, + launch.faceVolumes, + launch.geometryTotals, + launch.faceBendingEnergies, + launch.faceMeanCurvatures, + launch.faceNormals, + launch.occurrenceForces, + launch.sampleSurfaceMeasures, + launch.sampleMeanCurvatures, + launch.sampleNormals, + launch.sampleBendingEnergies, + launch.statusDiagnostics}) + clear(handle); + const std::size_t sampleCount = + static_cast(launch.evaluatedFaceCount) * + kQuadratureSampleCount; + for (std::size_t sample = 0; sample < sampleCount; ++sample) + { + write_double(launch.sampleSurfaceMeasures, sample, 1.0); + write_double(launch.sampleNormals, sample * 3U + 2U, 1.0); + } + for (std::size_t face = 0; + face < static_cast(launch.faceCount); ++face) + write_double(launch.faceNormals, face * 3U + 2U, 1.0); + switch (membraneCorruption) + { + case MembraneCorruption::None: + break; + case MembraneCorruption::Status: + write_i32(launch.statusDiagnostics, 0, + static_cast( + MembraneStatusCode::NonFiniteIntermediate)); + write_i32(launch.statusDiagnostics, 1, 0); + write_i32(launch.statusDiagnostics, 2, 1); + break; + case MembraneCorruption::NonFiniteOccurrence: + write_double(launch.occurrenceForces, 0, + std::numeric_limits::quiet_NaN()); + break; + case MembraneCorruption::ZeroSurfaceMeasure: + write_double(launch.sampleSurfaceMeasures, 0, 0.0); + break; + case MembraneCorruption::NonUnitNormal: + write_double(launch.sampleNormals, 2, 2.0); + break; + case MembraneCorruption::MismatchedTotals: + write_double(launch.geometryTotals, 0, 1.0); + break; + } + return DriverStatus{}; + }; ops.synchronize = [this]() { ++synchronizeCalls; if (synchronizeCalls == failSynchronizeCall) @@ -654,6 +741,105 @@ TEST(CudaMeshStateCoreTest, } } +TEST(CudaMeshStateCoreTest, + CandidateMembraneReturnsUnscatteredDiagnosticsAndCanRollback) +{ + FakeDevice device; + const RegularMeshPack pack = make_geometry_pack(); + auto created = create_mesh_state_core(device.operations(), pack); + ASSERT_NE(created.state, nullptr); + auto state = CudaMeshStateFactory::create( + std::move(created.state), created.report, + []() { return DriverStatus{}; }); + ASSERT_TRUE(state->prepare_candidate(pack.acceptedCoordinates, 2).ok()); + const MembraneCandidateResult result = + state->compute_candidate_membrane(); + ASSERT_TRUE(result.ok()) << result.error.message; + EXPECT_EQ(result.faceAreas.size(), pack.faceCount); + EXPECT_EQ(result.faceBendingEnergies.size(), pack.faceCount); + EXPECT_EQ(result.faceNormals.size(), pack.faceCount * 3U); + EXPECT_EQ(result.occurrenceForces.size(), + pack.evaluatedFaceCount * kRegularControlCount * 9U); + EXPECT_EQ(result.sampleSurfaceMeasures.size(), + pack.evaluatedFaceCount * kQuadratureSampleCount); + EXPECT_EQ(state->report().phase, TransactionPhase::Validated); + EXPECT_TRUE(state->rollback().ok()); + EXPECT_EQ(state->report().phase, TransactionPhase::IdleAccepted); +} + +TEST(CudaMeshStateCoreTest, + CandidateMembraneFailuresAndMalformedOutputsPreserveAcceptedState) +{ + const std::vector corruptions{ + FakeDevice::MembraneCorruption::Status, + FakeDevice::MembraneCorruption::NonFiniteOccurrence, + FakeDevice::MembraneCorruption::ZeroSurfaceMeasure, + FakeDevice::MembraneCorruption::NonUnitNormal, + FakeDevice::MembraneCorruption::MismatchedTotals, + }; + for (const FakeDevice::MembraneCorruption corruption : corruptions) + { + SCOPED_TRACE(static_cast(corruption)); + FakeDevice device; + const RegularMeshPack pack = make_geometry_pack(); + auto created = create_mesh_state_core(device.operations(), pack); + ASSERT_NE(created.state, nullptr); + const DeviceBufferHandle acceptedHandle = + created.state->accepted_coordinate_handle_for_testing(); + const std::vector acceptedBytes = + device.memory.at(acceptedHandle); + auto state = CudaMeshStateFactory::create( + std::move(created.state), created.report, + []() { return DriverStatus{}; }); + ASSERT_TRUE(state->prepare_candidate(pack.acceptedCoordinates, 2).ok()); + device.membraneCorruption = corruption; + const MembraneCandidateResult failed = + state->compute_candidate_membrane(); + EXPECT_EQ(failed.error.code, DeviceStateErrorCode::CandidateFailed); + if (corruption == FakeDevice::MembraneCorruption::Status) + { + EXPECT_EQ(failed.status, + MembraneStatusCode::NonFiniteIntermediate); + EXPECT_EQ(failed.failedSample, 1u); + } + EXPECT_EQ(state->report().phase, TransactionPhase::Failed); + EXPECT_EQ(state->report().lastOutcome, TransactionOutcome::Failed); + EXPECT_EQ(device.memory.at(acceptedHandle), acceptedBytes); + ASSERT_TRUE(state->recover().ok()); + EXPECT_EQ(state->report().phase, TransactionPhase::IdleAccepted); + EXPECT_EQ(device.memory.at(acceptedHandle), acceptedBytes); + } + + FakeDevice device; + const RegularMeshPack pack = make_geometry_pack(); + auto created = create_mesh_state_core(device.operations(), pack); + ASSERT_NE(created.state, nullptr); + auto state = CudaMeshStateFactory::create( + std::move(created.state), created.report, + []() { return DriverStatus{}; }); + ASSERT_TRUE(state->prepare_candidate(pack.acceptedCoordinates, 2).ok()); + device.failMembraneCall = device.membraneCalls + 1; + MembraneCandidateResult failed = state->compute_candidate_membrane(); + EXPECT_EQ(failed.error.code, DeviceStateErrorCode::CandidateFailed); + ASSERT_TRUE(state->recover().ok()); + + ASSERT_TRUE(state->prepare_candidate(pack.acceptedCoordinates, 2).ok()); + device.failMembraneCall = 0; + device.failCopyToHostCall = device.copyToHostCalls + 1; + failed = state->compute_candidate_membrane(); + EXPECT_EQ(failed.error.code, DeviceStateErrorCode::TransferFailed); + ASSERT_TRUE(state->recover().ok()); + + ASSERT_TRUE(state->prepare_candidate(pack.acceptedCoordinates, 2).ok()); + device.failCopyToHostCall = 0; + device.failSynchronizeCall = device.synchronizeCalls + 1; + failed = state->compute_candidate_membrane(); + EXPECT_EQ(failed.error.code, + DeviceStateErrorCode::SynchronizationFailed); + EXPECT_EQ(state->report().phase, TransactionPhase::Failed); + ASSERT_TRUE(state->recover().ok()); +} + TEST(CudaMeshStateCoreTest, UpdateAllocationFailurePreservesResidentState) { FakeDevice device; diff --git a/tests/test_cuda_mesh_state_inventory.py b/tests/test_cuda_mesh_state_inventory.py index 19d9e2c..2d631d0 100644 --- a/tests/test_cuda_mesh_state_inventory.py +++ b/tests/test_cuda_mesh_state_inventory.py @@ -23,13 +23,15 @@ def test_public_contract_and_core_have_required_anchors(self): "rollback()", "GeometryCandidateResult", "compute_candidate_geometry", + "MembraneCandidateResult", + "compute_candidate_membrane", ): self.assertIn(anchor, public) self.assertIn("DeviceOperations", (ROOT / "include/cuda/detail/Cuda_mesh_state_core.hpp").read_text()) self.assertIn("MemoryBudgetExceeded", core) self.assertIn("topology replacement requires fresh dependent generations", core) - def test_step_is_geometry_only_and_not_production_routed(self): + def test_step_contains_membrane_formula_but_is_not_production_routed(self): paths = [ ROOT / "include/cuda/Cuda_mesh_state.hpp", ROOT / "include/cuda/detail/Cuda_regular_geometry_cpu.hpp", @@ -40,9 +42,12 @@ def test_step_is_geometry_only_and_not_production_routed(self): text = "\n".join(path.read_text() for path in paths) self.assertIn("regular_geometry_kernel", text) self.assertIn("deterministic_geometry_reduction_kernel", text) + self.assertIn("regular_membrane_kernel", text) + self.assertIn("deterministic_membrane_reduction_kernel", text) + self.assertIn("occurrenceForces", text) self.assertNotIn("Compute_energy_and_force", text) self.assertNotIn("#include \"mesh/", text) - self.assertNotIn("force scatter", text.lower()) + self.assertNotIn("forceTotal", text) def test_make_targets_are_explicit_and_mutually_exclusive(self): makefile = (ROOT / "Makefile").read_text() @@ -156,6 +161,55 @@ def test_runner_rejects_false_green_geometry(self): missing["geometry_cases"].pop("curved") self.assertFalse(module.geometry_complete(missing)) + def test_runner_rejects_false_green_membrane(self): + path = ROOT / "scripts/run_cuda_mesh_state_report.py" + spec = importlib.util.spec_from_file_location("mesh_state_runner_membrane", path) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + complete_case = { + "pass": True, + "cpu_parity": True, + "repeatable": True, + "structured_degeneracy": True, + "recoverable": True, + "permutation_equal": True, + "max_abs_error": 1.0e-12, + } + complete = { + "membrane_repeatable": True, + "membrane_degeneracy_handled": True, + "membrane_max_abs_error": 1.0e-12, + "membrane_cases": { + name: complete_case.copy() + for name in module.REQUIRED_MEMBRANE_CASES + }, + } + self.assertTrue(module.membrane_complete(complete)) + for key, invalid in ( + ("membrane_repeatable", False), + ("membrane_degeneracy_handled", False), + ("membrane_max_abs_error", 1.0001e-12), + ("membrane_max_abs_error", "0"), + ): + with self.subTest(key=key, invalid=invalid): + report = complete.copy() + report[key] = invalid + self.assertFalse(module.membrane_complete(report)) + for name in module.REQUIRED_MEMBRANE_CASES: + for key, invalid in ( + ("pass", False), + ("cpu_parity", False), + ("repeatable", False), + ("structured_degeneracy", False), + ("recoverable", False), + ("permutation_equal", False), + ("max_abs_error", 1.0001e-12), + ): + with self.subTest(case=name, key=key): + report = json.loads(json.dumps(complete)) + report["membrane_cases"][name][key] = invalid + self.assertFalse(module.membrane_complete(report)) + def test_committed_rtx_evidence_meets_exit_gate(self): evidence = json.loads( (ROOT / "analysis/cuda_mesh_state_report_rtx4050.json").read_text() @@ -166,6 +220,9 @@ def test_committed_rtx_evidence_meets_exit_gate(self): self.assertEqual(evidence["iterations"], 20) self.assertLessEqual(evidence["geometry_max_abs_error"], 1.0e-12) self.assertTrue(evidence["geometry_repeatable"]) + self.assertLessEqual(evidence["membrane_max_abs_error"], 1.0e-12) + self.assertTrue(evidence["membrane_repeatable"]) + self.assertTrue(evidence["membrane_degeneracy_handled"]) self.assertEqual( set(evidence["geometry_cases"]), { @@ -185,6 +242,25 @@ def test_committed_rtx_evidence_meets_exit_gate(self): self.assertTrue(case["ghost_zero"]) self.assertTrue(case["degenerate_zero"]) self.assertTrue(case["permutation_equal"]) + self.assertEqual( + set(evidence["membrane_cases"]), + { + "natural", + "permuted", + "curved", + "boundary_ghost", + "degenerate", + "production_cpu", + }, + ) + for case in evidence["membrane_cases"].values(): + self.assertTrue(case["pass"]) + self.assertTrue(case["cpu_parity"]) + self.assertTrue(case["repeatable"]) + self.assertTrue(case["structured_degeneracy"]) + self.assertTrue(case["recoverable"]) + self.assertTrue(case["permutation_equal"]) + self.assertLessEqual(case["max_abs_error"], 1.0e-12) self.assertTrue(evidence["no_warm_allocations"]) self.assertTrue(evidence["transfers_complete"]) self.assertTrue(evidence["closed"]) diff --git a/tests/test_cuda_regular_membrane_cpu.cpp b/tests/test_cuda_regular_membrane_cpu.cpp new file mode 100644 index 0000000..91c66c1 --- /dev/null +++ b/tests/test_cuda_regular_membrane_cpu.cpp @@ -0,0 +1,127 @@ +#include "cuda/Cuda_mesh_pack.hpp" +#include "cuda/detail/Cuda_regular_membrane_cpu.hpp" +#include "mesh/Mesh.hpp" + +#include + +#include +#include + +namespace +{ +using namespace slimed::cuda_residency; +using namespace slimed::cuda_residency::detail; + +TEST(CudaRegularMembraneCpuTest, + PackedOracleMatchesProductionRegularFormulaAtEveryOutputLevel) +{ + Param param; + param.VERBOSE_MODE = false; + param.boundaryCondition = BoundaryType::Periodic; + param.sideX = 40.0; + param.sideY = 10.0 * std::sqrt(3.0) / 2.0 * param.lFace; + param.kCurv = 2.75; + param.uSurf = 1.25; + param.uVol = 0.85; + Mesh mesh(param); + ::testing::internal::CaptureStdout(); + mesh.setup_flat(); + ::testing::internal::GetCapturedStdout(); + for (Vertex &vertex : mesh.vertices) + { + const double index = static_cast(vertex.index); + vertex.coord.set(2, 0, 0.015 * std::sin(0.31 * index) + + 0.006 * std::cos(0.73 * index)); + vertex.coordPrev = vertex.coord; + vertex.coordRef = vertex.coord; + } + mesh.calculate_element_area_volume(); + + RegularMeshPackRequest request; + request.generations = {1, 1, 1, 1, 1}; + const RegularMeshPackResult packed = + build_regular_mesh_pack(mesh, request); + ASSERT_TRUE(packed.ok()) << packed.error.message; + const RegularMembraneCpuResult oracle = evaluate_regular_membrane_cpu( + packed.pack, packed.pack.acceptedCoordinates); + ASSERT_TRUE(oracle.ok()) << oracle.message; + + ASSERT_EQ(oracle.occurrenceForces.size(), + packed.pack.evaluatedFaceCount * kRegularControlCount * 9U); + ASSERT_EQ(oracle.sampleSurfaceMeasures.size(), + packed.pack.evaluatedFaceCount * kQuadratureSampleCount); + for (std::size_t evaluated = 0; + evaluated < static_cast(packed.pack.evaluatedFaceCount); + ++evaluated) + { + Face &face = mesh.faces[static_cast( + packed.pack.evaluatedFaceIds[evaluated])]; + std::vector coordinates; + coordinates.reserve(face.oneRingVertices.size()); + for (const int source : face.oneRingVertices) + coordinates.push_back(mesh.vertices[source].coord); + double meanCurvature = 0.0; + double bendingEnergy = 0.0; + Matrix normal = mat_calloc(3, 1); + Matrix bending = mat_calloc(kRegularControlCount, 3); + Matrix area = mat_calloc(kRegularControlCount, 3); + Matrix volume = mat_calloc(kRegularControlCount, 3); + mesh.element_energy_force_regular( + coordinates, face, face.spontCurvature, meanCurvature, normal, + bendingEnergy, bending, area, volume, true); + + EXPECT_NEAR(oracle.faceBendingEnergies[face.index], bendingEnergy, + 1.0e-11); + EXPECT_NEAR(oracle.faceMeanCurvatures[face.index], meanCurvature, + 1.0e-11); + for (std::size_t axis = 0; axis < 3; ++axis) + EXPECT_NEAR(oracle.faceNormals[face.index * 3 + axis], + normal.get(static_cast(axis), 0), 1.0e-11); + for (std::size_t local = 0; local < kRegularControlCount; ++local) + { + const std::size_t base = + (evaluated * kRegularControlCount + local) * 9U; + for (std::size_t axis = 0; axis < 3; ++axis) + { + EXPECT_NEAR(oracle.occurrenceForces[base + axis], + bending.get(static_cast(local), + static_cast(axis)), + 1.0e-10); + EXPECT_NEAR(oracle.occurrenceForces[base + 3 + axis], + area.get(static_cast(local), + static_cast(axis)), + 1.0e-10); + EXPECT_NEAR(oracle.occurrenceForces[base + 6 + axis], + volume.get(static_cast(local), + static_cast(axis)), + 1.0e-10); + } + } + } +} + +TEST(CudaRegularMembraneCpuTest, DegenerateSamplesReturnStructuredStatus) +{ + RegularMeshPack pack; + pack.vertexCount = kRegularControlCount; + pack.faceCount = 1; + pack.evaluatedFaceCount = 1; + pack.evaluatedFaceIds = {0}; + pack.oneRingSourceIds.resize(kRegularControlCount); + for (std::size_t local = 0; local < kRegularControlCount; ++local) + pack.oneRingSourceIds[local] = static_cast(local); + pack.evaluatedFaceSpontaneousCurvature = {0.0}; + pack.quadratureCoefficients.assign(kQuadratureSampleCount, 1.0 / 3.0); + pack.shapeWeights.assign(kQuadratureSampleCount * kShapeRowCount * + kRegularControlCount, + 0.0); + const std::vector coordinates(kRegularControlCount * 3U, 2.0); + const RegularMembraneCpuResult result = + evaluate_regular_membrane_cpu(pack, coordinates); + EXPECT_FALSE(result.ok()); + EXPECT_EQ(result.status, RegularMembraneStatus::DegenerateSample); + EXPECT_EQ(result.failedEvaluatedFace, 0u); + EXPECT_EQ(result.failedSample, 0u); +} + +} // namespace