update commits in hdf5_loop branch from v2.0 #5
Conversation
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
The nikkhoo_walter.mod file is a compiled artifact that should not be tracked in version control. It was accidentally added when checking out MATLAB reference files from the origin/SEAS_BP5_nikkhoo branch.
Our Fortran code calls tdstress_hs (half-space), but the previous MATLAB test used TDstressFS (full-space), which explains value mismatch. This test will confirm if the expected values (7.064e-4, 2.113e-4) come from TDstressHS half-space implementation. Need to run: matlab -nodisplay -r "test_points_8_9_HS; exit"
Modified TDstressHS.m to output individual contributions: - Main Dislocation (full-space triangular dislocation) - Harmonic Function (correction term) - Image Dislocation (mirror for half-space) - Total (sum of all three) This allows direct comparison with Fortran output to identify which component(s) differ and cause NaN in Points 8 & 9. Files: - TDstressHS.m: Added fprintf statements for each contribution - test_contributions_HS.m: Test script (nested functions prevent use) - COMPARISON_INSTRUCTIONS.md: Detailed instructions for running comparison Expected outcome: Identify whether MATLAB avoids singularities in Main/Image dislocation or gets finite result through cancellation.
Fixed critical bug where intermediate variables were declared but never calculated, causing NaN for ALL points including the triangle center. Variables that were declared but uninitialized: - W2 = W * W - Wr = W * r - W2r = W2 * r - Wr3 = W * r3 - W2r2 = W2 * r2 - rz = r * z - r2z2 = r2 * z2 - r3z = r3 * z These calculations match the MATLAB code (TDstressHS.m lines 614-618). The bug was introduced when singularity handling was removed - these essential calculations were accidentally deleted along with the checks. Result: Using uninitialized variables in formulas caused NaN for every point, not just edge cases. This explains why even the center of the triangle returned NaN. Added test files: - test_center.f90: Test specifically for triangle center - test_multiple_points.f90: Test suite for various point locations - test_center.sh: Build script Fix verified by matching MATLAB calculation structure.
CRITICAL BUG: Array indices were incorrect when translating from MATLAB's 2D arrays to Fortran's 3D arrays. MATLAB uses 2D points: p = [p(1), p(2)] = [y_coord, z_coord] Fortran uses 3D points: p = [p(1), p(2), p(3)] = [x_coord, y_coord, z_coord] Mapping: - MATLAB's p(1) → Fortran's p(2) (y-coordinate) - MATLAB's p(2) → Fortran's p(3) (z-coordinate) BEFORE (WRONG): denominator = (p2(2)-p3(2))*(p1(2)-p3(2)) + (p3(2)-p2(2))*(p1(3)-p3(3)) a = ((p2(2)-p3(2))*(x-p3(2)) + (p3(2)-p2(2))*(y-p3(3))) / denominator b = ((p3(2)-p1(2))*(x-p3(2)) + (p1(2)-p3(2))*(y-p3(3))) / denominator AFTER (CORRECT, matching MATLAB): denominator = (p2(3)-p3(3))*(p1(2)-p3(2)) + (p3(2)-p2(2))*(p1(3)-p3(3)) a = ((p2(3)-p3(3))*(x-p3(2)) + (p3(2)-p2(2))*(y-p3(3))) / denominator b = ((p3(3)-p1(3))*(x-p3(2)) + (p1(2)-p3(2))*(y-p3(3))) / denominator This bug caused incorrect barycentric coordinates, classifying the triangle center as b=0.0 (on edge) instead of b≈0.333 (inside triangle), leading to trimode=0 and NaN results via casez_log. MATLAB reference: TDstressHS.m lines 457-460
Fixed issue where Points 4, 5, 12, and 15 were returning NaN after the barycentric coordinate fix. PROBLEM: The Fortran code had overly complex edge detection that checked if barycentric coordinates were within [0,1] bounds with tolerance: if (abs(a) < TOL .and. b >= -TOL .and. b <= 1+TOL .and. c >= -TOL .and. c <= 1+TOL) This caused points that should be trimode=1 (first configuration) to be incorrectly classified as trimode=0 (edge case), leading to NaN results. SOLUTION: Simplified to match MATLAB's logic exactly (TDstressHS.m:467-470): MATLAB: trimode(a==0 & b>=0 & c>=0) = 0 Fortran: if (abs(a) < TOL .and. b >= 0 .and. c >= 0) then trimode = 0 Key changes: - Removed [0,1] bounds checking on b and c - Only check if a≈0 and b,c are non-negative (matching MATLAB) - Same for the other two edge cases This allows points with barycentric coordinates like (1.5, 0.0, -0.5) to be correctly classified as trimode=-1 (outside) instead of incorrectly triggering trimode=0 (edge). Affected points: - Point 4: (7.0, -1.0, -5.0) - Point 5: (-7.0, -1.0, -5.0) - Point 12: (1.0, -1.0, -1.0) - Point 15: (1.0, -1.0, -8.0) Added test_problem_points.f90 to specifically test these cases.
Points 4, 5, 12, and 15 return NaN after the three bug fixes. This appears to be CORRECT behavior - these points are geometrically singular (on or near triangle edges in the triangle plane). Files added: 1. test_matlab_points_4_5_12_15.m - MATLAB verification script Tests the same points with MATLAB TDstressHS to confirm behavior 2. debug_specific_points.f90 - Fortran debug program Analyzes why these points are singular 3. SINGULAR_POINTS_EXPLANATION.md - Comprehensive documentation Explains: - Why these points are singular (geometric analysis) - Coordinate system transformations (EFCS → TDCS) - The z≠0 override logic and when it applies - Physical interpretation of singularities - How to verify with MATLAB - What to do if finite values are needed Point analysis: - Point 4: (7.0, -1.0, -5.0) → x_td≈0, in triangle plane, on edge - Point 5: (-7.0, -1.0, -5.0) → x_td≈0, in triangle plane, on edge - Point 12: (1.0, -1.0, -1.0) → projects to edge/vertex - Point 15: (1.0, -1.0, -8.0) → projects to edge/vertex These points have trimode=0 (edge case) with z=0 (x_td=0, in plane), so the z≠0 override doesn't apply. They remain classified as singular edges, which triggers casez_log=TRUE → NaN. Next step: Run MATLAB verification to confirm MATLAB also returns NaN.
Found major bug in coord_trans and matrix construction that caused incorrect TDCS coordinates and wrong trimode classification. BUG: Matrix orientation mismatch between coord_trans and MATLAB ==================================================================== MATLAB uses different matrix orientations for different operations: 1. CoordTrans (coordinate transformation): - A = [Vnorm Vstrike Vdip]' (transposed, vectors as ROWS) - r = A * vector (matrix-vector multiplication) 2. TensTrans (tensor transformation): - Uses [Vnorm Vstrike Vdip] (NOT transposed, vectors as COLUMNS) BEFORE (WRONG): - Fortran A had vectors as COLUMNS: A(:,1)=vnorm, A(:,2)=vstrike, A(:,3)=vdip - coord_trans used: matmul(transpose(A), vector) = matmul(A_rows, vector) - Result: A_columns -> transpose -> A_rows (correct) BUT inconsistent with tdsetup_s AFTER (CORRECT): - Fortran A has vectors as ROWS: A(1,:)=vnorm, A(2,:)=vstrike, A(3,:)=vdip - coord_trans uses: matmul(A, vector) with A already having rows - tens_trans uses: transpose(A) to get columns for tensor transformation - Matches MATLAB exactly Changes: 1. sub_nikkhoo.f90:182-184 - A matrix construction (columns -> rows) 2. sub_nikkhoo.f90:443 - coord_trans (remove transpose, use A directly) 3. sub_nikkhoo.f90:281 - tens_trans call (add transpose(A)) 4. sub_nikkhoo.f90:505-507 - angsetup_fsc_s A matrix (columns -> rows) IMPACT: This bug caused completely wrong TDCS coordinates, leading to: - Wrong barycentric coordinates - Wrong trimode classification - Points 4,5,12,15 incorrectly classified as singular (trimode=0) - NaN results for valid points that should be finite Expected results after fix: - Point 4: Exx = 0.000829157341339727 - Point 5: Exx = 0.00114439668841158 - Point 12: Exx = 0.00441202690885827 - Point 15: Exx = -0.000914111766849476 Added verify_trimode.f90 to test these points.
Created complete Python package translating MATLAB TDstressFS.m and TDstressHS.m for triangular dislocation stress/strain calculations. Files added: ============ Core modules: - td_utils.py: Helper functions for coordinate/tensor transformations, trimodefinder - ang_dislocation.py: Angular dislocation strain calculations - tdstress_fs.py: Full-space triangular dislocation solution - tdstress_hs.py: Half-space solution (main + image dislocation) Package infrastructure: - __init__.py: Package initialization, exports main functions - setup.py: Installation script for pip install - requirements.txt: Dependencies (numpy) - README.md: Comprehensive documentation with examples - test_tdstress.py: Test script with reference values Features: ========= 1. Full-Space Solution (tdstress_fs): - Calculate stress/strain for triangular dislocations in infinite elastic medium - Handles three configurations based on point location - Returns NaN for singular points (on edges in plane) 2. Half-Space Solution (tdstress_hs): - Calculates main dislocation + image dislocation - Accounts for free surface effects - Note: Harmonic function contribution simplified (placeholder) 3. Helper Functions: - coord_trans: Coordinate system transformations (EFCS ↔ TDCS) - tens_trans: Tensor component transformations - trimodefinder: Determine configuration (inside/outside/edge) - ang_dis_strain: Angular dislocation strain - td_setup_s: Setup and transform for angular dislocations 4. Well-documented: - Docstrings for all functions - README with examples and usage - Test script with reference values Translation notes: ================== - Direct translation from MATLAB to Python/NumPy - Maintains same coordinate systems (EFCS, TDCS, ADCS) - Uses numpy arrays instead of MATLAB matrices - Function signatures match MATLAB where possible - Handles scalar and array inputs Known limitations: ================== 1. Harmonic function (AngSetupFSC_S) simplified - returns zeros Full implementation would require: - Angular dislocation pairs on each triangle side - Free surface correction calculations - Complex transformation logic 2. Singular points return NaN (correct but may need regularization) 3. No special optimizations for large arrays yet Reference: ========== Nikkhoo M. and Walter T.R., 2015. Triangular dislocation: An analytical, artefact-free solution. Geophysical Journal International. Original MATLAB: Mehdi Nikkhoo ([email protected])
Implemented all components for the free surface correction in half-space: New modules: - ang_dislocation_fsc.py: Free surface correction strains for angular dislocation (AngDisStrainFSC with all 6 strain components) - ang_setup_fsc.py: Setup and calculation for angular dislocation pairs on each triangle side (AngSetupFSC_S) Updates: - tdstress_hs.py: Complete implementation of tdstress_harfunc with harmonic function contribution - __init__.py: Updated exports and version to 2.0.0 - README.md: Updated documentation to reflect complete implementation - test_tdstress.py: Fixed imports to work with package structure Implementation details: - Translated ~600 lines of complex MATLAB formulas for 6 strain components - Handles two configurations based on (beta*y1A) >= 0 - Calculates contributions from all three triangle sides (P1P2, P2P3, P3P1) - Transforms between EFCS and ADCS coordinate systems - Computes both stress and strain tensors Note: Additional testing and validation against MATLAB reference values is needed to verify numerical accuracy. Refs: Nikkhoo & Walter (2015) - Triangular dislocation method
Updated test_tdstress.py to test all 15 points from the Fortran test suite: - Expanded test cases from 5 to 15 points - Added proper handling for points without reference values - Updated half-space test to use all 5 reference points - Improved output formatting for better readability Created TEST_RESULTS.md documenting complete test results: - Full-space results for all 15 points - Half-space results for 5 key points - Detailed comparison showing Python has same bugs as original Fortran - Identified Bug #2 (barycentric coordinates) and Bug #3 (edge detection) - Documented spurious NaNs at points 4, 5, 12, 15 - Recommended fixes matching Fortran bug patches Test results confirm: - Point 1 (center): Wrong value due to barycentric bug - Points 4, 5, 12, 15: Spurious NaNs (should be finite) - Point 11 (on edge): Correct NaN (true singularity) - Point 2 (surface): Numerical overflow needs handling Python implementation needs same bug fixes as Fortran code.
Corrected all 15 test point coordinates to match the Fortran arrays: - Point 1: (-1/3, -1/3, -3.0) - Point 2 (center): (-1/3, -1/3, -14/3) - Point 3: (-1/3, -1/3, -6.0) - Points 4-15: Fixed to match exact Fortran x[i], y[i], z[i] arrays Previous coordinates were incorrect and didn't match the Fortran test. Updated TEST_RESULTS.md with corrected results: - 6 points now return spurious NaNs (points 4, 5, 6, 9, 12, 15) - All are well away from triangle edges/vertices - Confirms same bugs as original Fortran code Test results show Python implementation has identical bugs to unpatched Fortran: Bug #2 (barycentric) and Bug #3 (edge detection)
Updated test_tdstress.py with: - All 15 test points matching Fortran test exactly - Complete expected Exx values from corrected Fortran reference - Comprehensive error reporting for both TDstressFS and TDstressHS Updated TEST_RESULTS.md with: - Full results table for all 15 points (both full-space and half-space) - Detailed analysis showing 0 PASS, 9 FAIL, 6 NaN - Critical issues: spurious NaNs at 6 valid points, large errors everywhere - Root cause analysis linking to same bugs as original Fortran (Bug #2, #3, #4) - Specific fix recommendations for each bug - Runtime warnings documentation Key findings: - Python exhibits SAME BUGS as original unpatched Fortran code - Bug #2 (barycentric coordinates): center point 88% error - Bug #3 (edge detection): 6 spurious NaNs at valid distant points - Bug #4 (matrix orientation): likely causing large errors (up to 3961%) Ready for bug fixes to be applied following Fortran corrections.
Analysis comparing MATLAB and Python code reveals Bug #4: - Python incorrectly uses A.T (transpose) in 4 locations - MATLAB uses A directly (no transpose) in same locations - Affects TDstress_HarFunc and AngSetupFSC_S Critical issues found: 1. tdstress_hs.py:163 - slip vector transformation uses A.T (should be A) 2. ang_setup_fsc.py:92 - point transformation uses A.T (should be A) 3. ang_setup_fsc.py:95 - side vector transformation uses A.T (should be A) 4. ang_setup_fsc.py:103 - slip vector transformation uses A.T (should be A) Impact: Coordinate transformations reversed, causing errors up to 3961% Root cause: Misunderstanding of MATLAB matrix construction conventions Documented in MATRIX_TRANSPOSE_ISSUES.md with: - Side-by-side MATLAB vs Python comparison - Explanation of correct usage in TDstressFS - Root cause analysis - Fix recommendations for all 4 locations
Fixed 4 locations where Python incorrectly used A.T instead of A: 1. tdstress_hs.py:163 - slip vector transformation (TDCS to EFCS) 2. ang_setup_fsc.py:92 - point transformation (EFCS to ADCS) 3. ang_setup_fsc.py:95 - side vector transformation (EFCS to ADCS) 4. ang_setup_fsc.py:103 - slip vector transformation (EFCS to ADCS) Changes: - Removed .T from coord_trans() calls to match MATLAB reference - MATLAB uses A directly when A has unit vectors as columns - Python was incorrectly transposing these transformations Impact on test results: - Slight improvements in some half-space calculations - Point 1: 32.2% → 29.9% error (improved) - Point 2: 87.0% → 86.7% error (slight improvement) - Point 3: 3955% → 3954% error (marginal) Note: Large errors still persist, indicating additional bugs remain (Bug #2: barycentric coordinates, Bug #3: edge detection)
Created comprehensive comparison showing: - Before vs after results for all 15 test points - Improvements: Points 1, 2, 11, 13 (marginal, 2-10% reduction) - Degradations: Points 7, 8 (got worse) - No change: Points 3, 10, 14 and all NaN points Key findings: - Transpose fixes necessary but insufficient - Errors still massive (30-3954%) - Point 2 (center): still 86.7% error (Bug #2) - 6 spurious NaNs persist (Bug #3) - Point 3: 3954% error indicates additional formula bugs Conclusion: Bug #2 (barycentric) and Bug #3 (edge detection) must be fixed next to make meaningful progress.
Ignore: - __pycache__/ and compiled Python files - Virtual environments - IDE files (.vscode, .idea) - OS files (.DS_Store) - Fortran compiled files
Critical bug found in td_utils.py:144-146: - Python extracts p1[:2] = [x_TDCS, y_TDCS] (WRONG) - MATLAB uses p1(2:3) = [y_TDCS, z_TDCS] (CORRECT) Impact: - Barycentric coordinates calculated in wrong 2D plane (x-y instead of y-z) - All trimode classifications incorrect - Center point has 88% error due to wrong barycentric coords Root cause: - MATLAB: passes p1(2:3), p2(2:3), p3(2:3) (last 2 elements) - Python: extracts p1[:2], p2[:2], p3[:2] (first 2 elements) Fix needed: Change p1[:2] to p1[1:3] in three locations This is Bug #2 from Fortran analysis (wrong barycentric indices). Documented in BUG2_ARRAY_EXTRACTION.md with detailed analysis.
Audited all major array/matrix extractions between MATLAB and Python: ✅ CORRECT (4 implementations): 1. SideVec(1:2) extraction in AngSetupFSC_S 2. Matrix A construction in TDSetupS 3. Matrix B construction in TDSetupS 4. Barycentric formula structure ❌ INCORRECT (1 implementation - Bug #2): - Vertex array extraction in trimodefinder - MATLAB: passes p1(2:3) = [y_TDCS, z_TDCS] - Python: extracts p1[:2] = [x_TDCS, y_TDCS] - Projects onto WRONG plane (x-y instead of y-z) - Causes 88% error at triangle center Fix: Change p1[:2] to p1[1:3] in td_utils.py:144-146 This confirms Bug #2 is the array indexing error, matching the original Fortran bug diagnosis.
Created new MatlabScript folder with comprehensive stiffness matrix computation: compute_trigreen_matrix.m: - Computes elastic Green's function (stiffness matrix) for triangular meshes - Uses Nikkhoo & Walter (2015) TDstressHS method - Outputs binary files compatible with Fortran TriGreen code - Supports parallel distribution across multiple processors - Configurable material properties (mu, nu) and slip components Features: - Reads GTS mesh format - Distributes work across ncore processors - Outputs trigreen_<id>.bin files (Nt × Ncell values per file) - Creates position.bin with element centroids - Automatic load balancing for uneven cell distribution example_run.m: - Demonstrates 6 usage examples - Shows default and custom material properties - Examples for different slip configurations - Output directory specification - Single and multi-processor configurations README.md: - Complete documentation with usage instructions - File format specifications (input GTS, output binary) - Performance considerations and memory requirements - Comparison with Fortran implementation - Troubleshooting guide and verification methods test_mesh.gts: - Simple 8-element test mesh for immediate testing - 3×3 grid of vertices forming 8 triangular elements Output format matches Fortran calc_trigreen.f90: - trigreen_<id>.bin: Binary files with Nt*Ncell double values - position.bin: Element centroid positions - Units: Bar (0.1 MPa)
Fixed critical array extraction bug in td_utils.py:144-146: - Changed p1[:2] to p1[1:3] to extract [y_TDCS, z_TDCS] - Changed p2[:2] to p2[1:3] - Changed p3[:2] to p3[1:3] Added debug scripts to verify barycentric calculation: - debug_barycentric.py: Shows transformation and barycentric coords - debug_formula.py: Manual formula trace showing correct result Verification: - Manual calculation with correct mapping gives a=b=c=0.333 ✓ - Formula is now correctly using [y_TDCS, z_TDCS] projection Note: Despite correct fix, test errors remain at 88%. This indicates additional bugs exist beyond the array extraction issue. The barycentric coordinates are now correct, but errors persist in other parts of the code. Next: Investigate why correct barycentric coords don't improve results. Likely issues in angular dislocation formulas or coordinate transformations.
Removed __pycache__/*.pyc files from git tracking: - These files were previously tracked before .gitignore was added - They are now properly ignored by .gitignore - Files remain locally but won't be tracked by git This prevents the 'uncommitted changes' warning from appearing when running tests that regenerate these cache files.
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
No description provided.