Fix singularity in tetrahedron evaluation. Disabled OpenMP in micromag

rasmusbj · rasmusbj · commit a1cfec2a121e · 2023-09-21T12:36:33.000+02:00
diff --git a/buildMagTenseMEX.m b/buildMagTenseMEX.m
@@ -115,7 +115,7 @@ function buildMagTenseMEX(USE_RELEASE,USE_CUDA,USE_CVODE)
 %% MagTenseLandauLifshitzSolver_mex
 Source_str = 'source/MagTenseMEX/MagTenseMEX/MagTenseLandauLifshitzSolver_mex.f90';
 mex_str = ['mex ' compiler_str ' ' Debug_flag ' ' MagTenseMicroMag_str ' ' DemagField_str ' ' TileDemagTensor_str ' ' NumericalIntegration_str ' ' CUDA_str ' ' CVODE_str ' ' MKL_str ' ' Source_str ' ' Options_str];
-eval(mex_str) 
+eval_MEX(mex_str)
 if ~USE_RELEASE
     pause(pause_time)
     movefile(['MagTenseLandauLifshitzSolver_mex.mex' MEX_str '64.pdb'],['matlab/MEX_files/MagTenseLandauLifshitzSolver_mex.mex' MEX_str '64.pdb']);
@@ -127,7 +127,7 @@ function buildMagTenseMEX(USE_RELEASE,USE_CUDA,USE_CVODE)
 MagTenseMicroMag_str        = ['-Lsource/MagTenseMicroMag/' build_str_NO_CUDA_MagTenseMicroMag '/ ' MagTenseMicroMag_lib_str ' -Isource/MagTenseMicroMag/' build_str_MagTenseMicroMag '/'];
 Source_str = 'source/MagTenseMEX/MagTenseMEX/MagTenseLandauLifshitzSolver_mex.f90';
 mex_str = ['mex ' compiler_str ' ' Debug_flag ' ' MagTenseMicroMag_str ' ' DemagField_str ' ' TileDemagTensor_str ' ' NumericalIntegration_str ' ' CVODE_str ' ' MKL_str ' ' Source_str ' ' Options_str];
-eval(mex_str) 
+eval_MEX(mex_str)
 if ~USE_RELEASE
     pause(pause_time)
     movefile(['MagTenseLandauLifshitzSolver_mex.mex' MEX_str '64.pdb'],['matlab/MEX_files/MagTenseLandauLifshitzSolverNoCUDA_mex.mex' MEX_str '64.pdb']);
@@ -139,7 +139,7 @@ function buildMagTenseMEX(USE_RELEASE,USE_CUDA,USE_CVODE)
 %% IterateMagnetization_mex
 Source_str = 'source/MagTenseMEX/MagTenseMEX/IterateMagnetization_mex.f90';
 mex_str = ['mex ' compiler_str ' ' Debug_flag ' ' DemagField_str ' ' TileDemagTensor_str ' ' NumericalIntegration_str ' ' Source_str ' ' Options_str];
-eval(mex_str) 
+eval_MEX(mex_str)
 if ~USE_RELEASE
     pause(pause_time)
     movefile(['IterateMagnetization_mex.mex' MEX_str '64.pdb'],['matlab/MEX_files/IterateMagnetization_mex.mex' MEX_str '64.pdb']);
@@ -150,7 +150,7 @@ function buildMagTenseMEX(USE_RELEASE,USE_CUDA,USE_CVODE)
 %% getHFromTiles_mex
 Source_str = 'source/MagTenseMEX/MagTenseMEX/getHFromTiles_mex.f90';
 mex_str = ['mex ' compiler_str ' ' Debug_flag ' ' DemagField_str ' ' TileDemagTensor_str ' ' NumericalIntegration_str ' ' Source_str ' ' Options_str];
-eval(mex_str) 
+eval_MEX(mex_str)
 if ~USE_RELEASE
     pause(pause_time)
     movefile(['getHFromTiles_mex.mex' MEX_str '64.pdb'],['matlab/MEX_files/getHFromTiles_mex.mex' MEX_str '64.pdb']);
@@ -161,7 +161,7 @@ function buildMagTenseMEX(USE_RELEASE,USE_CUDA,USE_CVODE)
 %% getNFromTile_mex
 Source_str = 'source/MagTenseMEX/MagTenseMEX/getNFromTile_mex.f90';
 mex_str = ['mex ' compiler_str ' ' Debug_flag ' ' DemagField_str ' ' NumericalIntegration_str ' ' TileDemagTensor_str ' ' Source_str ' ' Options_str];
-eval(mex_str) 
+eval_MEX(mex_str)
 if ~USE_RELEASE
     pause(pause_time)
     movefile(['getNFromTile_mex.mex' MEX_str '64.pdb'],['matlab/MEX_files/getNFromTile_mex.mex' MEX_str '64.pdb']);
@@ -172,7 +172,7 @@ function buildMagTenseMEX(USE_RELEASE,USE_CUDA,USE_CVODE)
 %% getMagForce_mex
 Source_str = 'source/MagTenseMEX/MagTenseMEX/getMagForce_mex.f90';
 mex_str = ['mex ' compiler_str ' ' Debug_flag ' ' MagneticForceIntegrator_str ' ' DemagField_str ' ' NumericalIntegration_str ' ' TileDemagTensor_str ' ' Source_str ' ' Options_str];
-eval(mex_str) 
+eval_MEX(mex_str)
 if ~USE_RELEASE
     pause(pause_time)
     movefile(['getMagForce_mex.mex' MEX_str '64.pdb'],['matlab/MEX_files/getMagForce_mex.mex' MEX_str '64.pdb']);
@@ -181,3 +181,26 @@ function buildMagTenseMEX(USE_RELEASE,USE_CUDA,USE_CVODE)
 movefile(['getMagForce_mex.mex' MEX_str '64'],['matlab/MEX_files/getMagForce_mex.mex' MEX_str '64']);
     
 end
+
+function eval_MEX(mex_str)
+
+try
+    eval(mex_str);
+catch ME
+    if (strcmp(ME.message(91:117),'mt : general error c101008d'))
+        fail_mex = true; 
+        while fail_mex
+            try 
+                disp('Microsoft manifest tool error - retrying')
+                eval(mex_str); 
+                fail_mex = false; 
+            catch
+                continue
+            end
+        end
+    else
+        disp(ME.message)
+        rethrow(ME)
+    end
+end
+end
diff --git a/source/DemagField/DemagField/IterateMagnetSolution.f90 b/source/DemagField/DemagField/IterateMagnetSolution.f90
@@ -204,6 +204,12 @@ subroutine iterateMagnetization( tiles, n, stateFunction, n_stf, T, err_max, max
                         
                         pts = pts + tiles(i)%offset     !< Finally add the offset to translate the circ piece
                             
+                    case(tileTypeTetrahedron)
+                        !! No rotation is needed as the center of the tetrahedron is with respect to the vertices
+                        pts(1) = (tiles(i)%vert(1,1) + tiles(i)%vert(1,2) + tiles(i)%vert(1,3) + tiles(i)%vert(1,4))/4.0
+                        pts(2) = (tiles(i)%vert(2,1) + tiles(i)%vert(2,2) + tiles(i)%vert(2,3) + tiles(i)%vert(2,4))/4.0
+                        pts(3) = (tiles(i)%vert(3,1) + tiles(i)%vert(3,2) + tiles(i)%vert(3,3) + tiles(i)%vert(3,4))/4.0
+
             
                     case default
                 
diff --git a/source/MagTenseMicroMag/LandauLifshitzEquationSolver.f90 b/source/MagTenseMicroMag/LandauLifshitzEquationSolver.f90
@@ -699,7 +699,7 @@ subroutine updateDemagfield( problem, solution)
             beta = 0.0
             !Hmx = Kxx * Mx
             call gemv( problem%Kxx, solution%Mx_s, solution%HmX, alpha, beta )
-            
+                                   
             beta = 1.0
             !Hmx = Hmx + Kxy * My
             call gemv( problem%Kxy, solution%My_s, solution%HmX, alpha, beta )
@@ -914,8 +914,9 @@ subroutine ComputeDemagfieldTensor( problem )
 
         CALL SYSTEM_CLOCK(c1)
  
-        call mkl_set_num_threads(problem%nThreadsMatlab)
-        call omp_set_num_threads(problem%nThreadsMatlab)
+        !call mkl_set_num_threads(problem%nThreadsMatlab)
+        !call omp_set_num_threads(problem%nThreadsMatlab)
+        call omp_set_num_threads(1)
                
         if ( problem%grid%gridType .eq. gridTypeUniform ) then
             
@@ -984,10 +985,10 @@ subroutine ComputeDemagfieldTensor( problem )
                 call displayGUIMessage( 'Averaging the N_tensor not supported for this tile type' )
             endif
     
-            !$OMP PARALLEL shared(problem)
-            !$omp do private(ind, tile, H, Nout)
-            
+            !$OMP PARALLEL shared(problem) private(ind, indx_ele, tile, H, Nout, pts_arr, i)
+                        
             !for each element find the tensor for all evaluation points (i.e. all elements)
+            !$omp do
             do i=1,nx
                 !Setup template tile
                 tile(1)%tileType = 5 !(for tetrahedron)
@@ -1001,7 +1002,13 @@ subroutine ComputeDemagfieldTensor( problem )
                 allocate(Nout(1,ntot,3,3))
                 allocate(H(ntot,3))
                 
-                call getFieldFromTiles( tile, H, problem%grid%pts, 1, ntot, Nout, .false. )
+                allocate(pts_arr(ntot,3))
+                pts_arr(:,1) =  problem%grid%pts(:,1)
+                pts_arr(:,2) =  problem%grid%pts(:,2)
+                pts_arr(:,3) =  problem%grid%pts(:,3)
+                
+                !call getFieldFromTiles( tile, H, problem%grid%pts, 1, ntot, Nout, .false. )
+                call getFieldFromTiles( tile, H, pts_arr, 1, ntot, Nout, .false. )
                     
                 !Copy Nout into the proper structure used by the micro mag model
                 ind = i
@@ -1022,13 +1029,15 @@ subroutine ComputeDemagfieldTensor( problem )
                 problem%Kzz(:,ind) = sngl(Nout(1,:,3,3))
                 
                 !Clean up
+                deallocate(pts_arr)
                 deallocate(Nout)
                 deallocate(H)
             enddo
-            
             !$omp end do
+            
             !$OMP END PARALLEL
-        
+            
+            
         elseif ( problem%grid%gridType .eq. gridTypeUnstructuredPrisms ) then
             
             !$OMP PARALLEL shared(problem)
@@ -1129,8 +1138,7 @@ subroutine ComputeDemagfieldTensor( problem )
         
         !call displayGUIMessage( 'Kxx(1,1):' )
         !write (prog_str,'(f10.3)') problem%Kxx(1,1)
-        !call displayGUIMessage( prog_str )
-        
+        !call displayGUIMessage( prog_str )      
         
         !Write the demag tensors to disk if asked to do so            
         if ( problem%demagTensorReturnState .gt. DemagTensorReturnMemory ) then
diff --git a/source/MagTenseMicroMag/MicroMagParameters.f90 b/source/MagTenseMicroMag/MicroMagParameters.f90
@@ -105,7 +105,7 @@ module MicroMagParameters
         integer :: demag_approximation                          !> Flag for how to approximate the demagnetization tensor as specified in the parameters below
         integer :: demagTensorReturnState                       !> Flag describing how or if the demag tensor should be returned
         integer :: demagTensorLoadState                         !> Flag describing how or if to load the demag tensor (from disk e.g.)
-        integer*4 :: nThreadsMatlab                             !> Number of threads to use in the OpenMP demag tensor allocation
+        integer :: nThreadsMatlab                             !> Number of threads to use in the OpenMP demag tensor allocation
         integer,dimension(3) :: N_ave                           !> Number of points to average the demag tensor in in the recieving tile, N_ave(1) = N_x etc
         character*256 :: demagTensorFileOut, demagTensorFileIn  !> Filename (including path) for output (input) of demag tensor if it is to be returned as a file (demagTensorReturnState >2 and the value is equal to the length of the file including path)
         
diff --git a/source/TileDemagTensor/TileDemagTensor/TileTriangle.f90 b/source/TileDemagTensor/TileDemagTensor/TileTriangle.f90
@@ -1,5 +1,6 @@
 module TileTriangle
     use SPECIALFUNCTIONS
+    use, intrinsic :: ieee_arithmetic
     implicit none
     
     real,parameter :: pi=3.14159265359
@@ -169,8 +170,8 @@ function Nxz( r, l, h )
     
     ! Returns the Nxz tensor component in the local coordinate system
         Nxz = -1./(4.*pi) * ( F_Nxz(r,h,l,h) - F_Nxz(r,0.,l,h) - ( G_Nxz(r,h) - G_Nxz(r,0.) ) )
-        
-    end function Nxz
+    
+        end function Nxz
 
     function F_Nxz( r, yp, l, h )
     real :: F_Nxz
@@ -189,6 +190,10 @@ function G_Nxz( r, yp )
     real,dimension(3),intent(in) :: r
     real,intent(in) :: yp
         G_Nxz = atanh( ( r(2)-yp )/ (sqrt(r(1)**2+(r(2)-yp)**2+r(3)**2)) )
+        
+        if (.not.ieee_is_finite(G_Nxz)) then     ! variable is not finite
+            G_Nxz = 0;
+        endif
     end function G_Nxz
 
     !---------------------------------------------------------------------------    
@@ -223,6 +228,10 @@ function L_Nyz( r, xp )
     real,dimension(3),intent(in) :: r
     real,intent(in) :: xp
         L_Nyz = atanh( (r(1) - xp) / (sqrt((r(1)-xp)**2+r(2)**2+r(3)**2)) )
+        
+        if (.not.ieee_is_finite(L_Nyz)) then     ! variable is not finite
+            L_Nyz = 0;
+        endif
     end function L_Nyz
 
     !---------------------------------------------------------------------------    
