Merge lp:~r-nelson/fluidity/parallel_sigma into lp:fluidity

[Rhodri Nelson]

Fixed a bug with Sigma coordinates in parallel. Instead of using a pre-existing chain as a dummy chain to create the initial mesh, a specific dummy chain is now created. This ensures consistency across all processors.

Also, updated a couple of interfaces to be more in line with the fluidity style.

[Simon Funke (simon-funke)]

This merge request has empty? What happened?

review: Needs Information

[Lawrence Mitchell (wence)]

I think the relevant differences in the branch are already in trunk (see r3545)

[Adam Candy (asc)]

That was my conclusion too.
It'd be good to check that there isn't some extra code that was meant to be included in this merge - Rhodri? Let us know when you're back online - さっぽろはどうですか？　おげんきで！

[Rhodri Nelson (r-nelson)]

Something strange seems to have happened, which was probably my mistake. I modified my trunk and then made a branch from it but for some reason none of the changes were ported into the branch. Maybe it's because I forgot to commit them locally or something? I can't access Rossby at the moment though so can't get hold of the code. I'll update the other merge request later though as I have that on my laptop.

[lp:~r-nelson/fluidity/parallel_sigma]

3551. By Rhodri Nelson on 2011-09-08

Updating some local changes.

[Rhodri Nelson (r-nelson)]

The relevant changes have now been uploaded to Launchpad.

Unmerged revisions

3551. By Rhodri Nelson on 2011-09-08

Updating some local changes.

=== modified file 'horizontal_adaptivity/Combine_Meshes.F90'
--- horizontal_adaptivity/Combine_Meshes.F90	2011-08-01 16:03:12 +0000
+++ horizontal_adaptivity/Combine_Meshes.F90	2011-09-08 00:21:25 +0000
@@ -22,27 +22,32 @@
   contains
 
   subroutine combine_z_meshes(h_mesh, z_meshes, out_mesh, &
-    full_shape, mesh_name, option_path, sl)
-  !! Given the h_mesh and a z_mesh under each node of it combines these
-  !! into a full horiz+vertic. mesh
-  type(vector_field), intent(inout):: h_mesh
-  type(vector_field), dimension(:), intent(in):: z_meshes
-  type(vector_field), intent(out):: out_mesh
-  type(element_type), intent(in):: full_shape
-  !! the name of the topol. mesh to be created, not the coordinate field
-  character(len=*), intent(in):: mesh_name, option_path
-  logical, intent(in), optional :: sl
-  logical :: sigma_layers
+    full_shape, mesh_name, option_path, sigma_layers)
+
+    !! Given the h_mesh and a z_mesh under each node of it combines these
+    !! into a full horiz+vertic. mesh
+    type(vector_field), intent(inout):: h_mesh
+    type(vector_field), dimension(:), intent(in):: z_meshes
+    type(vector_field), intent(out):: out_mesh
+    type(element_type), intent(in):: full_shape
+    !! the name of the topol. mesh to be created, not the coordinate field
+    character(len=*), intent(in):: mesh_name, option_path
+    logical, intent(in), optional :: sigma_layers
+
+    ! The sigma_layers logical input is optional so create another logical for use
+    ! in this routine  
+    logical :: sl
+    type(vector_field) :: dummy_column
   
     type(csr_sparsity):: out_columns
     type(mesh_type):: mesh
     integer, dimension(:), allocatable:: no_hanging_nodes
     integer:: column, total_out_nodes, total_out_elements, z_elements, last_seen
 
-    if (present(sl)) then
-      sigma_layers=sl
+    if (present(sigma_layers)) then
+      sl=sigma_layers
     else
-      sigma_layers=.false.
+      sl=.false.
     end if
     
     allocate(no_hanging_nodes(1:node_count(h_mesh)))
@@ -95,14 +100,22 @@
     out_mesh%option_path=""
     out_mesh%mesh%periodic = mesh_periodic(h_mesh)
 
+    ! If we have sigma layers we want to create a dummy chain in order to create a mesh
+    ! with the desired topological properties. The nodal coordinates will then later
+    ! be over written with the correct ones. Note that all chains have the same amount of
+    ! layers so we will just use the first chain to create the dummy mesh.
+    if (sl) then
+      call create_dummy_column(z_meshes(1), dummy_column)
+    end if
+
     last_seen = 0
     do column=1,node_count(h_mesh)
       if (node_owned(h_mesh, column)) then
-        if (sigma_layers) then
+        if (sl) then
           ! If we have sigma layers we will first use one chain to chain to create a dummy
           ! 'mesh' that has the correct topological properties. We will later over write the
           !  vector field with the correct one. We always have chain '1', so we'll use that now.
-          call append_to_structures(column, z_meshes(1), h_mesh, out_mesh, last_seen)
+          call append_to_structures(column, dummy_column, h_mesh, out_mesh, last_seen)
         else
           call append_to_structures(column, z_meshes(column), h_mesh, out_mesh, last_seen)
         end if
@@ -115,6 +128,7 @@
         last_seen = last_seen + no_hanging_nodes(column)+1
       end if
     end do
+    if (sl) call deallocate(dummy_column)
     assert(all(out_mesh%mesh%columns>0))
       
     call create_columns_sparsity(out_columns, out_mesh%mesh)
@@ -130,13 +144,23 @@
 
     ! If we have sigma layers we now populate the vector field with the actual
     ! node positions.
-    if (sigma_layers) then
+    if (sl) then
       last_seen = 0
       do column=1,node_count(h_mesh)
         if (node_owned(h_mesh, column)) then
           call append_to_structures(column, z_meshes(column), h_mesh, out_mesh, last_seen)
+        else
+          ! for non-owned columns we reserve node numbers, 
+          ! but don't fill in out_mesh positions yet
+          ! note that in this way out_mesh will obtain the same halo ordering
+          ! convention as h_mesh
+          out_mesh%mesh%columns(last_seen+1:last_seen+no_hanging_nodes(column)+1) = column
+          last_seen = last_seen + no_hanging_nodes(column)+1
         end if
       end do
+      if (associated(h_mesh%mesh%halos)) then
+        call halo_update(out_mesh)
+      end if
     end if
 
     call deallocate(out_columns)
@@ -145,27 +169,32 @@
   
 
   subroutine combine_r_meshes(shell_mesh, r_meshes, out_mesh, &
-    full_shape, mesh_name, option_path, sl)
-  !! Given the shell_mesh and a r_mesh under each node of it combines these
-  !! into a full layered mesh
-  type(vector_field), intent(inout):: shell_mesh
-  type(vector_field), dimension(:), intent(in):: r_meshes
-  type(vector_field), intent(out):: out_mesh
-  type(element_type), intent(in):: full_shape
-  !! the name of the topol. mesh to be created, not the coordinate field
-  character(len=*), intent(in):: mesh_name, option_path
-  logical, intent(in), optional :: sl
-  logical :: sigma_layers
+    full_shape, mesh_name, option_path, sigma_layers)
+
+    !! Given the shell_mesh and a r_mesh under each node of it combines these
+    !! into a full layered mesh
+    type(vector_field), intent(inout):: shell_mesh
+    type(vector_field), dimension(:), intent(in):: r_meshes
+    type(vector_field), intent(out):: out_mesh
+    type(element_type), intent(in):: full_shape
+    !! the name of the topol. mesh to be created, not the coordinate field
+    character(len=*), intent(in):: mesh_name, option_path
+    logical, intent(in), optional :: sigma_layers
+
+    ! The sigma_layers logical input is optional so create another logical for use
+    ! in this routine  
+    logical :: sl
+    type(vector_field) :: dummy_column
   
     type(csr_sparsity):: out_columns
     type(mesh_type):: mesh
     integer, dimension(:), allocatable:: no_hanging_nodes
     integer:: column, total_out_nodes, total_out_elements, r_elements, last_seen
 
-    if (present(sl)) then
-      sigma_layers=sl
+    if (present(sigma_layers)) then
+      sl=sigma_layers
     else
-      sigma_layers=.false.
+      sl=.false.
     end if
     
     allocate(no_hanging_nodes(1:node_count(shell_mesh)))
@@ -217,15 +246,23 @@
     out_mesh%mesh%option_path=option_path
     out_mesh%option_path=""
     out_mesh%mesh%periodic = mesh_periodic(shell_mesh) ! Leave this here for now
+
+    ! If we have sigma layers we want to create a dummy chain in order to create a mesh
+    ! with the desired topological properties. The nodal coordinates will then later
+    ! be over written with the correct ones. Note that all chains have the same amount of
+    ! layers so we will just use the first chain to create the dummy mesh.
+    if (sl) then
+      call create_dummy_column(r_meshes(1), dummy_column)
+    end if
     
     last_seen = 0
     do column=1,node_count(shell_mesh)
       if (node_owned(shell_mesh, column)) then
-        if (sigma_layers) then
+        if (sl) then
           ! If we have sigma layers we will first use one chain to chain to create a dummy
           ! 'mesh' that has the correct topological properties. We will later over write the
           !  vector field with the correct one. We always have chain '1', so we'll use that now.
-          call append_to_structures_radial(column, r_meshes(1), out_mesh, last_seen)
+          call append_to_structures_radial(column, r_meshes(column), out_mesh, last_seen, dummy_column)
         else
           call append_to_structures_radial(column, r_meshes(column), out_mesh, last_seen)
         end if
@@ -238,6 +275,7 @@
         last_seen = last_seen + no_hanging_nodes(column)+1
       end if
     end do
+    if (sl) call deallocate(dummy_column)
     assert(all(out_mesh%mesh%columns>0))
       
     call create_columns_sparsity(out_columns, out_mesh%mesh)
@@ -251,16 +289,25 @@
 
     call generate_layered_mesh(out_mesh, shell_mesh)
 
-    if (sigma_layers) then
+    ! If we have sigma layers we now populate the vector field with the actual
+    ! node positions.
+    if (sl) then
       last_seen = 0
       do column=1,node_count(shell_mesh)
         if (node_owned(shell_mesh, column)) then
-          ! If we have sigma layers we will first use one chain to chain to create a dummy
-          ! 'mesh' that has the correct topological properties. We will later over write the
-          !  vector field with the correct one. We always have chain '1', so we'll use that now.
           call append_to_structures_radial(column, r_meshes(column), out_mesh, last_seen)
+        else
+          ! for non-owned columns we reserve node numbers, 
+          ! but don't fill in out_mesh positions yet
+          ! note that in this way out_mesh will obtain the same halo ordering
+          ! convention as h_mesh
+          out_mesh%mesh%columns(last_seen+1:last_seen+no_hanging_nodes(column)+1) = column
+          last_seen = last_seen + no_hanging_nodes(column)+1
         end if
       end do
+      if (associated(shell_mesh%mesh%halos)) then
+        call halo_update(out_mesh)
+      end if
     end if
 
     call deallocate(out_columns)
@@ -291,25 +338,55 @@
     
   end subroutine append_to_structures
 
-  subroutine append_to_structures_radial(column, r_mesh, out_mesh, last_seen)
+  subroutine append_to_structures_radial(column, r_mesh, out_mesh, last_seen, dummy_column)
     integer, intent(in) :: column
     type(vector_field), intent(in) :: r_mesh
     type(vector_field), intent(inout) :: out_mesh
     integer, intent(inout) :: last_seen
+    type(vector_field), intent(in), optional :: dummy_column
 
     integer :: j
     integer :: r_dim
 
     real, dimension(mesh_dim(out_mesh)) :: pos
 
+    logical :: sigma_layers
+    real, dimension(3) :: r
+    real :: mod_r, theta, phi, height
+
+    if (present(dummy_column)) then
+      sigma_layers=.true.
+    else
+      sigma_layers=.false.
+    end if
+
     r_dim = mesh_dim(out_mesh)
 
-    do j=1,node_count(r_mesh)
-      last_seen = last_seen + 1
-      out_mesh%mesh%columns(last_seen)=column
-      pos(1:r_dim) = node_val(r_mesh, j)             ! note that in z-dir extrusion case only the z value was loaded into the out mesh from the z mesh. Here however,
-      call set(out_mesh, last_seen, pos)
-    end do
+    if (sigma_layers) then
+      assert(mesh_dim(out_mesh)==3)
+      r=node_val(r_mesh, 1)
+      mod_r=sqrt(dot_product(r,r))
+      
+      theta=acos(r(3)/mod_r)
+      phi=atan2(r(2),r(1))
+ 
+      do j=1,node_count(r_mesh)
+        last_seen = last_seen + 1
+        out_mesh%mesh%columns(last_seen)=column
+        height=node_val(dummy_column,1,j)
+        pos(1)=(mod_r+height)*sin(theta)*cos(phi)
+        pos(2)=(mod_r+height)*sin(theta)*sin(phi)
+        pos(3)=(mod_r+height)*cos(theta)
+        call set(out_mesh, last_seen, pos)
+      end do
+    else
+      do j=1,node_count(r_mesh)
+        last_seen = last_seen + 1
+        out_mesh%mesh%columns(last_seen)=column
+        pos(1:r_dim) = node_val(r_mesh, j)      ! note that in z-dir extrusion case only the z value was loaded into the out mesh from the z mesh. Here however,
+        call set(out_mesh, last_seen, pos)
+      end do
+    end if
     
   end subroutine append_to_structures_radial
 
@@ -405,5 +482,48 @@
     
   end subroutine derive_extruded_l2_node_halo
 
+  subroutine create_dummy_column(z_mesh, dummy_column)
+  ! Subroutine to create a dummy column when using sigma layers
+
+    type(vector_field), intent(in) :: z_mesh
+    type(vector_field), intent(inout) :: dummy_column
+
+    integer :: node
+
+    type(mesh_type) :: mesh
+    type(element_type) :: oned_shape
+    type(quadrature_type) :: oned_quad
+    integer :: quadrature_degree
+    integer :: ele
+    integer, parameter :: loc=2
+
+    integer :: elements
+    type(rlist):: depths
+
+    call get_option("/geometry/quadrature/degree", quadrature_degree)
+    oned_quad = make_quadrature(vertices=loc, dim=1, degree=quadrature_degree)
+    oned_shape = make_element_shape(vertices=loc, dim=1, degree=1, quad=oned_quad)
+    call deallocate(oned_quad)
+  
+    do node=1,node_count(z_mesh)
+      call insert(depths, float(node-1)*10.0)
+    end do
+    elements=depths%length-1
+
+    call allocate(mesh, elements+1, elements, oned_shape, "DummyMesh")
+    call deallocate(oned_shape)
+    do ele=1,elements
+      mesh%ndglno((ele-1) * loc + 1: ele*loc) = (/ele, ele+1/)
+    end do
+
+    call allocate(dummy_column, 1, mesh, "DummyMeshCoordinates")
+    call deallocate(mesh)
+
+    call set(dummy_column, 1, (/0.0/))
+    do node=1, elements+1
+      call set(dummy_column, node,  (/ pop(depths) /))
+    end do
+
+  end subroutine create_dummy_column
 
 end module hadapt_combine_meshes

=== modified file 'horizontal_adaptivity/Extrude.F90'
--- horizontal_adaptivity/Extrude.F90	2011-08-01 16:03:12 +0000
+++ horizontal_adaptivity/Extrude.F90	2011-09-08 00:21:25 +0000
@@ -155,7 +155,7 @@
 
     ! combine the 1d vertical meshes into a full mesh
     call combine_z_meshes(h_mesh, z_meshes, out_mesh, &
-       full_shape, mesh_name, option_path, sigma_layers)
+       full_shape, mesh_name, option_path, sigma_layers=sigma_layers)
        
     do column=1, node_count(h_mesh)
       if (.not. node_owned(h_mesh, column)) cycle

=== modified file 'horizontal_adaptivity/Extrude_radially.F90'
--- horizontal_adaptivity/Extrude_radially.F90	2011-08-01 16:03:12 +0000
+++ horizontal_adaptivity/Extrude_radially.F90	2011-09-08 00:21:25 +0000
@@ -285,7 +285,7 @@
 
     ! combine the 1d radial meshes into a full mesh
     call combine_r_meshes(shell_mesh, r_meshes, out_mesh, &
-       full_shape, mesh_name, option_path, sigma_layers)
+       full_shape, mesh_name, option_path, sigma_layers=sigma_layers)
 
 !     vtk_write_state(filename, index, model, state, write_region_ids, stat)
 !     call vtk_write_fields("extruded_mesh", 0, out_mesh, out_mesh%mesh, vfields=(/out_mesh/))