[cig-commits] [commit] master: in case of fluid core, sort all fluid elements by radius (0487a1a)

[cig_noreply at geodynamics.org]

Repository : https://github.com/geodynamics/axisem

On branch  : master
Link       : https://github.com/geodynamics/axisem/compare/49cf2d34d650ddeadc90951f0a3884d9dc5aff2e...0487a1a65c9a4ff020c97672b15153d9ab7736fb

>---------------------------------------------------------------

commit 0487a1a65c9a4ff020c97672b15153d9ab7736fb
Author: martinvandriel <vandriel at erdw.ethz.ch>
Date:   Fri Sep 26 15:55:29 2014 +0200

    in case of fluid core, sort all fluid elements by radius
    
    to make sure that the domains are connected and every domain as maximum 8 neighbours


>---------------------------------------------------------------

0487a1a65c9a4ff020c97672b15153d9ab7736fb
 MESHER/parallelization.f90 | 82 ++++++++++++++++++++++++----------------------
 1 file changed, 43 insertions(+), 39 deletions(-)

diff --git a/MESHER/parallelization.f90 b/MESHER/parallelization.f90
index 7722768..ff9165d 100644
--- a/MESHER/parallelization.f90
+++ b/MESHER/parallelization.f90
@@ -451,7 +451,7 @@ subroutine domain_decomposition_theta_r(attributed, nprocl, nthetal, nrl, &
   integer, intent(in)       :: nprocl, nthetal, nrl
   integer, intent(in)       :: nelmax, nelmax_fluid, nelmax_solid
   
-  integer                   :: itheta, iitheta, iel
+  integer                   :: itheta, iitheta, iel, nel_fluid_theta
   integer                   :: irad, iproc
   integer                   :: mycount, nicb, ncmb
   integer                   :: iprocb(2), mycountb(2), j1, j2
@@ -510,44 +510,6 @@ subroutine domain_decomposition_theta_r(attributed, nprocl, nthetal, nrl, &
   ! **************** END OF INNER CUBE****************
 
 
-  ! sort inner core elements according to radius
-  ! Using same stupid choice as in inner core decomposition 
-
-  if (solid_domain(ndisc)) then 
-     do itheta = 0, nthetal-1
-        allocate(inner_core_buf(central_count(itheta)))
-        allocate(inner_core_r(central_count(itheta)))
-        inner_core_buf(:) = thetaslel_solid(1:central_count(itheta),itheta)
-
-        do iel = 1, central_count(itheta)
-           inner_core_r(iel) = rcom(inner_core_buf(iel))
-        enddo
-  
-        call mergesort_3(inner_core_r, il=inner_core_buf, p=4)
-        
-        thetaslel_solid(1:central_count(itheta),itheta) = inner_core_buf(:)
-
-        deallocate(inner_core_buf)
-        deallocate(inner_core_r)
-     enddo
-  else
-     do itheta = 0, nthetal-1
-        allocate(inner_core_buf(central_count(itheta)))
-        allocate(inner_core_r(central_count(itheta)))
-        inner_core_buf(:) = thetaslel_fluid(1:central_count(itheta),itheta)
-
-        do iel = 1, central_count(itheta)
-           inner_core_r(iel) = rcom(inner_core_buf(iel))
-        enddo
-  
-        call mergesort_3(inner_core_r, il=inner_core_buf, p=4)
-        
-        thetaslel_fluid(1:central_count(itheta),itheta) = inner_core_buf(:)
-
-        deallocate(inner_core_buf)
-        deallocate(inner_core_r)
-     enddo
-  endif
   
   ! add the extra requirement that element iel to be in appropriate theta slice
   do itheta = 0, nthetal-1
@@ -657,6 +619,48 @@ subroutine domain_decomposition_theta_r(attributed, nprocl, nthetal, nrl, &
      stop
   endif
 
+
+  ! sort inner core elements according to radius
+  ! Using same stupid choice as in inner core decomposition 
+
+  if (solid_domain(ndisc)) then 
+     do itheta = 0, nthetal-1
+        allocate(inner_core_buf(central_count(itheta)))
+        allocate(inner_core_r(central_count(itheta)))
+        inner_core_buf(:) = thetaslel_solid(1:central_count(itheta),itheta)
+
+        do iel = 1, central_count(itheta)
+           inner_core_r(iel) = rcom(inner_core_buf(iel))
+        enddo
+  
+        call mergesort_3(inner_core_r, il=inner_core_buf, p=4)
+        
+        thetaslel_solid(1:central_count(itheta),itheta) = inner_core_buf(:)
+
+        deallocate(inner_core_buf)
+        deallocate(inner_core_r)
+     enddo
+  else
+     do itheta = 0, nthetal-1
+        nel_fluid_theta = sum(nel_fluid(itheta:itheta+nrl-1))
+        allocate(inner_core_buf(nel_fluid_theta))
+        allocate(inner_core_r(nel_fluid_theta))
+        inner_core_buf(:) = thetaslel_fluid(1:nel_fluid_theta,itheta)
+
+        do iel = 1, nel_fluid_theta
+           ! sort by radius, if radius is the same, theta makes the difference
+           inner_core_r(iel) = rcom(inner_core_buf(iel)) + 1e-10 * thetacom(inner_core_buf(iel))
+        enddo
+
+        call mergesort_3(inner_core_r, il=inner_core_buf, p=4)
+
+        thetaslel_fluid(1:nel_fluid_theta,itheta) = inner_core_buf(:)
+
+        deallocate(inner_core_buf)
+        deallocate(inner_core_r)
+     enddo
+  endif
+
   ! reset, as we have to touch each element again!
   attributed = .false.