← Back to branch summary

~fluidity-core/fluidity/tidal-velocity-bcs-doodson

« back to all changes in this revision

Viewing changes to assemble/Momentum_DG.F90

Merging in moving mesh checkpoint fix, so we can checkpoint our tidal runs.

Show diffs side-by-side

added added

removed removed

Lines of Context:
assemble/Momentum_DG.F90
47
47
  use boundary_conditions_from_options
48
48
  use solvers
49
49
  use dgtools
50
 
  use global_parameters, only: OPTION_PATH_LEN, FIELD_NAME_LEN
 
50
  use global_parameters, only: OPTION_PATH_LEN, FIELD_NAME_LEN, COLOURING_DG2, &
 
51
       COLOURING_DG0
51
52
  use coriolis_module
52
53
  use halos
53
54
  use sparsity_patterns
61
62
  use sparsity_patterns_meshes
62
63
  use colouring
63
64
  use Profiler
 
65
#ifdef _OPENMP
64
66
  use omp_lib
 
67
#endif
65
68
  use multiphase_module
66
69
 
67
70
 
220
223
    integer, dimension(:,:), allocatable :: velocity_bc_type
221
224
    integer, dimension(:), allocatable :: pressure_bc_type
222
225
    
223
 
 
224
226
    !! Sparsity for inverse mass
225
227
    type(csr_sparsity):: mass_sparsity
226
228
    
249
251
    real, dimension(u%dim) :: abs_wd_const
250
252
 
251
253
    !! 
252
 
    type(mesh_type) :: p0_mesh
253
 
    type(csr_sparsity), pointer :: dependency_sparsity
254
 
    type(scalar_field) :: node_colour
255
 
    type(integer_set), dimension(:), allocatable :: clr_sets
256
 
    integer :: clr, nnid, no_colours, len
257
 
    logical :: compact_stencil
 
254
    type(integer_set), dimension(:), pointer :: colours
 
255
    integer :: len, clr, nnid
 
256
    !! Is the transform_to_physical cache we prepopulated valid
 
257
    logical :: cache_valid
258
258
    integer :: num_threads
259
259
 
260
260
    ! Volume fraction fields for multi-phase flow simulation
617
617
      call remap_field(wettingdrying_alpha, alpha_u_field)
618
618
    end if
619
619
 
620
 
      call profiler_tic(u, "element_loop")
621
 
      
 
620
    call profiler_tic(u, "element_loop-omp_overhead")
622
621
 
623
622
#ifdef _OPENMP
624
 
      num_threads = omp_get_max_threads()
 
623
    num_threads = omp_get_max_threads()
625
624
#else 
626
 
      num_threads=1
 
625
    num_threads=1
627
626
#endif
628
627
 
629
 
    if(num_threads>1) then
630
 
 
631
 
      !! working out the options for different schemes of different terms
632
 
      call set_coriolis_parameters
633
 
      compact_stencil = have_option(trim(u%option_path)//&
634
 
            &"/prognostic/spatial_discretisation"//&
635
 
            &"/discontinuous_galerkin/viscosity_scheme"//&
636
 
            &"/interior_penalty") .or. &
637
 
            &have_option(trim(u%option_path)//&
638
 
            &"/prognostic/spatial_discretisation"//&
639
 
            &"/discontinuous_galerkin/viscosity_scheme"//&
640
 
            &"/compact_discontinuous_galerkin")
641
 
 
642
 
      compact_stencil=.false.
643
 
      !! generate the dual graph of the mesh
644
 
      p0_mesh = piecewise_constant_mesh(x%mesh, trim(x%name)//"P0Mesh")
645
 
 
646
 
       !! the sparse pattern of the dual graph.
647
 
       if (have_viscosity .and. (.not. compact_stencil)) then
648
 
          dependency_sparsity => get_csr_sparsity_secondorder(state, p0_mesh, p0_mesh)
649
 
       else
650
 
          dependency_sparsity =>get_csr_sparsity_compactdgdouble(state, p0_mesh)
651
 
       end if
652
 
       
653
 
       dependency_sparsity => get_csr_sparsity_secondorder(state, p0_mesh, p0_mesh)
654
 
 
655
 
       !! colouring dual graph according the sparsity pattern with greedy
656
 
       !! colouring algorithm
657
 
       !! "colours" is an array of type(integer_set)
658
 
       call colour_sparsity(dependency_sparsity, p0_mesh, node_colour, no_colours)
659
 
 
660
 
       if(.not. verify_colour_sparsity(dependency_sparsity, node_colour)) then
661
 
        FLAbort("The neighbours are using same colours, wrong!!.")
662
 
       endif 
663
 
 
664
 
       allocate(clr_sets(no_colours))
665
 
       clr_sets=colour_sets(dependency_sparsity, node_colour, no_colours)
666
 
       PRINT *, "colouring passed"
 
628
    if (have_viscosity) then
 
629
       call get_mesh_colouring(state, u%mesh, COLOURING_DG2, colours)
667
630
    else
668
 
       no_colours = 1
669
 
       allocate(clr_sets(no_colours))
670
 
       call allocate(clr_sets)
671
 
       do ELE=1,element_count(U)
672
 
          call insert(clr_sets(1), ele)
673
 
       end do
674
 
    end if
675
 
 
676
 
    colour_loop: do clr = 1, no_colours
677
 
      len = key_count(clr_sets(clr))
678
 
      !$OMP PARALLEL DO DEFAULT(SHARED) &
679
 
      !$OMP SCHEDULE(STATIC) &
680
 
      !$OMP PRIVATE(nnid, ele)
681
 
      element_loop: do nnid = 1, len 
682
 
       ele = fetch(clr_sets(clr), nnid)       
 
631
       call get_mesh_colouring(state, u%mesh, COLOURING_DG0, colours)
 
632
    end if
 
633
#ifdef _OPENMP
 
634
    cache_valid = prepopulate_transform_cache(X)
 
635
    assert(cache_valid)
 
636
    if (have_coriolis) then
 
637
       call set_coriolis_parameters
 
638
    end if
 
639
#endif
 
640
    call profiler_toc(u, "element_loop-omp_overhead")
 
641
    
 
642
    call profiler_tic(u, "element_loop")
 
643
 
 
644
    !$OMP PARALLEL DEFAULT(SHARED) &
 
645
    !$OMP PRIVATE(clr, nnid, ele, len)
 
646
 
 
647
    colour_loop: do clr = 1, size(colours) 
 
648
      len = key_count(colours(clr))
 
649
 
 
650
      !$OMP DO SCHEDULE(STATIC)
 
651
      element_loop: do nnid = 1, len
 
652
       ele = fetch(colours(clr), nnid)
683
653
       call construct_momentum_element_dg(ele, big_m, rhs, &
684
654
            & X, U, advecting_velocity, U_mesh, X_old, X_new, &
685
655
            & Source, Buoyancy, gravity, Abs, Viscosity, &
691
661
            & inverse_mass=inverse_mass, &
692
662
            & inverse_masslump=inverse_masslump, &
693
663
            & mass=mass, subcycle_m=subcycle_m)
694
 
      
695
664
      end do element_loop
696
 
      !$OMP END PARALLEL DO
 
665
      !$OMP END DO
697
666
 
698
667
    end do colour_loop
699
 
 
700
 
    call deallocate(clr_sets)
701
 
    deallocate(clr_sets)
702
 
 
703
 
    if(num_threads > 1) then
704
 
    call deallocate(node_colour)
705
 
    call deallocate(p0_mesh)
706
 
    endif
707
 
    
 
668
    !$OMP END PARALLEL
 
669
 
708
670
    call profiler_toc(u, "element_loop")
709
671
 
710
672
    if (have_wd_abs) then
893
855
    !Switch to select if we are assembling the primal or dual form
894
856
    logical :: primal
895
857
 
896
 
    ! In parallel, we only assemble the mass terms for the halo elements. All
897
 
    ! other terms are only assembled on elements we own.
898
 
    logical :: owned_element
 
858
    ! In parallel, we assemble terms on elements we own, and those in
 
859
    ! the L1 element halo
 
860
    logical :: assemble_element
899
861
 
900
862
    ! If on the sphere evaluate gravity direction at the gauss points
901
863
    logical :: on_sphere
941
903
    dg=continuity(U)<0
942
904
    p0=(element_degree(u,ele)==0)
943
905
    
944
 
    ! In parallel, we only construct the equations on elements we own.
945
 
    owned_element=element_owned(U,ele).or..not.dg
 
906
    ! In parallel, we construct terms on elements we own and those in
 
907
    ! the L1 element halo.
 
908
    assemble_element = .not.dg.or.element_neighbour_owned(U, ele)
946
909
 
947
910
    primal = .not.dg
948
911
    if(viscosity_scheme == CDG) primal = .true.
1050
1013
      deallocate(dnvfrac_t)
1051
1014
    end if
1052
1015
 
1053
 
    if ((have_viscosity).and.owned_element) then
 
1016
    if ((have_viscosity).and.assemble_element) then
1054
1017
      Viscosity_ele = ele_val(Viscosity,ele)
1055
1018
    end if
1056
1019
   
1057
 
    if (owned_element) then
 
1020
    if (assemble_element) then
1058
1021
       u_val = ele_val(u, ele)
1059
1022
    end if
1060
1023
 
1086
1049
       ! NOTE: this doesn't deal with mesh movement
1087
1050
       call addto(mass, u_ele, u_ele, Rho_mat)
1088
1051
    else
1089
 
      if(have_mass.and.owned_element) then
 
1052
      if(have_mass.and.assemble_element) then
1090
1053
        if(lump_mass) then        
1091
1054
          do dim = 1, u%dim
1092
1055
            big_m_diag_addto(dim, :loc) = big_m_diag_addto(dim, :loc) + l_masslump
1097
1060
          end do
1098
1061
        end if
1099
1062
      end if
1100
 
      if (move_mesh.and.owned_element) then
 
1063
      if (move_mesh.and.assemble_element) then
1101
1064
        ! In the unaccelerated form we solve:
1102
1065
        !  /
1103
1066
        !  |  N^{n+1} u^{n+1}/dt - N^{n} u^n/dt + ... = f
1122
1085
      end if
1123
1086
    end if
1124
1087
    
1125
 
    if(have_coriolis.and.(rhs%dim>1).and.owned_element) then
 
1088
    if(have_coriolis.and.(rhs%dim>1).and.assemble_element) then
1126
1089
      Coriolis_q=coriolis(ele_val_at_quad(X,ele))
1127
1090
    
1128
1091
      ! Element Coriolis parameter matrix.
1138
1101
      end if
1139
1102
    end if
1140
1103
 
1141
 
    if(have_advection.and.(.not.p0).and.owned_element) then
 
1104
    if(have_advection.and.(.not.p0).and.assemble_element) then
1142
1105
      ! Advecting velocity at quadrature points.
1143
1106
      U_nl_q=ele_val_at_quad(U_nl,ele)
1144
1107
 
1254
1217
 
1255
1218
    end if
1256
1219
 
1257
 
    if(have_source.and.acceleration.and.owned_element) then
 
1220
    if(have_source.and.acceleration.and.assemble_element) then
1258
1221
      ! Momentum source matrix.
1259
1222
      Source_mat = shape_shape(U_shape, ele_shape(Source,ele), detwei*Rho_q)
1260
1223
      if(lump_source) then
1271
1234
      end if
1272
1235
    end if
1273
1236
 
1274
 
    if(have_gravity.and.acceleration.and.owned_element) then
 
1237
    if(have_gravity.and.acceleration.and.assemble_element) then
1275
1238
      ! buoyancy
1276
1239
      if (on_sphere) then
1277
1240
      ! If were on a spherical Earth evaluate the direction of the gravity vector
1294
1257
      end if
1295
1258
    end if
1296
1259
 
1297
 
    if((have_absorption.or.have_vertical_stabilization.or.have_wd_abs) .and. (owned_element .or. pressure_corrected_absorption)) then
 
1260
    if((have_absorption.or.have_vertical_stabilization.or.have_wd_abs) .and. &
 
1261
         (assemble_element .or. pressure_corrected_absorption)) then
1298
1262
 
1299
1263
      absorption_gi=0.0
1300
1264
      tensor_absorption_gi=0.0
1381
1345
        if(lump_abs) then
1382
1346
 
1383
1347
          Abs_lump_sphere = sum(Abs_mat_sphere, 4)
1384
 
          if (owned_element) then          
 
1348
          if (assemble_element) then
1385
1349
            do dim = 1, U%dim
1386
1350
              do dim2 = 1, U%dim
1387
1351
                do i = 1, ele_loc(U, ele)
1415
1379
 
1416
1380
        else
1417
1381
 
1418
 
          if (owned_element) then
 
1382
          if (assemble_element) then
1419
1383
            do dim = 1, u%dim
1420
1384
              do dim2 = 1, u%dim
1421
1385
                big_m_tensor_addto(dim, dim2, :loc, :loc) = big_m_tensor_addto(dim, dim2, :loc, :loc) + &
1460
1424
        if(lump_abs) then        
1461
1425
          abs_lump = sum(Abs_mat, 3)
1462
1426
          do dim = 1, u%dim
1463
 
            if (owned_element) then
 
1427
            if (assemble_element) then
1464
1428
              big_m_diag_addto(dim, :loc) = big_m_diag_addto(dim, :loc) + dt*theta*abs_lump(dim,:)
1465
1429
              if(acceleration) then
1466
1430
                rhs_addto(dim, :loc) = rhs_addto(dim, :loc) - abs_lump(dim,:)*u_val(dim,:)
1481
1445
        else
1482
1446
      
1483
1447
          do dim = 1, u%dim
1484
 
            if (owned_element) then
 
1448
            if (assemble_element) then
1485
1449
              big_m_tensor_addto(dim, dim, :loc, :loc) = big_m_tensor_addto(dim, dim, :loc, :loc) + &
1486
1450
                & dt*theta*Abs_mat(dim,:,:)
1487
1451
              if(acceleration) then
1529
1493
    
1530
1494
    ! Viscosity.
1531
1495
    Viscosity_mat=0
1532
 
    if(have_viscosity.and.owned_element) then
 
1496
    if(have_viscosity.and.assemble_element) then
1533
1497
       if (primal) then
1534
1498
          do dim = 1, u%dim
1535
1499
             if(multiphase) then
1606
1570
       end if
1607
1571
    end if
1608
1572
 
1609
 
    if(have_surfacetension.and.(.not.p0).and.owned_element) then
 
1573
    if(have_surfacetension.and.(.not.p0).and.assemble_element) then
1610
1574
      if(integrate_surfacetension_by_parts) then
1611
1575
        tension = ele_val_at_quad(surfacetension, ele)
1612
1576
        
1624
1588
    ! Interface integrals
1625
1589
    !-------------------------------------------------------------------
1626
1590
    
1627
 
    if(dg.and.(have_viscosity.or.have_advection.or.have_pressure_bc).and.owned_element) then
 
1591
    if(dg.and.(have_viscosity.or.have_advection.or.have_pressure_bc).and.assemble_element) then
1628
1592
      neigh=>ele_neigh(U, ele)
1629
1593
      ! x_neigh/=t_neigh only on periodic boundaries.
1630
1594
      x_neigh=>ele_neigh(X, ele)
1654
1618
            ! Internal faces.
1655
1619
            face_2=ele_face(U, ele_2, ele)
1656
1620
        ! Check if face is turbine face (note: get_entire_boundary_condition only returns "applied" boundaries and we reset the apply status in each timestep)
1657
 
        elseif (velocity_bc_type(1,face)==4 .or. velocity_bc_type(1,face)==5) then  
1658
 
            face_2=face_neigh(turbine_conn_mesh, face)
1659
 
            turbine_face=.true.
 
1621
        elseif (velocity_bc_type(1,face)==4 .or. velocity_bc_type(1,face)==5) then
 
1622
           face_2=face_neigh(turbine_conn_mesh, face)
 
1623
           turbine_face=.true.
1660
1624
        else 
1661
1625
           ! External face.
1662
1626
           face_2=face
1793
1757
                        
1794
1758
                            ! Add BC into RHS
1795
1759
                            !
1796
 
                            rhs_addto(dim,:) = rhs_addto(dim,:) &
1797
 
                                & -matmul(Viscosity_mat(dim,:,start:finish), &
1798
 
                                & ele_val(velocity_bc,dim,face))
1799
 
                        
 
1760
                             rhs_addto(dim,:) = rhs_addto(dim,:) &
 
1761
                                  & -matmul(Viscosity_mat(dim,:,start:finish), &
 
1762
                                  & ele_val(velocity_bc,dim,face))
1800
1763
                            ! Ensure it is not used again.
1801
1764
                            Viscosity_mat(dim,:,start:finish)=0.0
1802
1765
                            
1837
1800
    ! Perform global assembly.
1838
1801
    !----------------------------------------------------------------------
1839
1802
    
1840
 
    if (owned_element) then
 
1803
    if (assemble_element) then
1841
1804
 
1842
1805
       ! add lumped terms to the diagonal of the matrix
1843
1806
       call add_diagonal_to_tensor(big_m_diag_addto, big_m_tensor_addto)
2103
2066
    if (boundary) then
2104
2067
       do dim=1,U%dim
2105
2068
          if (velocity_bc_type(dim,face)==1) then
2106
 
            dirichlet(dim)=.true.
 
2069
             dirichlet(dim)=.true.
2107
2070
          end if
2108
2071
       end do
2109
2072
       ! free surface b.c. is set for the 1st (normal) component
2235
2198
               ! Neumann boundary it's necessary to apply downwinding here
2236
2199
               ! to maintain the surface integral. Fortunately, since
2237
2200
               ! face_2==face for a boundary this is automagic.
2238
 
               
 
2201
 
2239
2202
               if (acceleration) then
2240
2203
                  rhs_addto(dim,u_face_l) = rhs_addto(dim,u_face_l) &
2241
2204
                       ! Outflow boundary integral.
2245
2208
               end if
2246
2209
               
2247
2210
            else
2248
 
               
 
2211
 
2249
2212
               rhs_addto(dim,u_face_l) = rhs_addto(dim,u_face_l) &
2250
2213
                    ! Outflow boundary integral.
2251
2214
                    -matmul(nnAdvection_out,face_val(U,dim,face))&
2253
2216
                    -matmul(nnAdvection_in,ele_val(velocity_bc,dim,face))
2254
2217
            end if
2255
2218
         end if
2256
 
       end do
 
2219
      end do
2257
2220
        
2258
2221
    end if
2259
2222
 
3135
3098
    
3136
3099
    ! we first work everything out for rows corresponding to the first component
3137
3100
    do ele=1, element_count(u)
3138
 
      ! we only have to provide nnz for owned rows
3139
 
      ! eventhough we do non-local assembly. The owner
 
3101
      ! we only have to provide nnz for owned rows. The owner
3140
3102
      ! therefore needs to specify the correct nnzs including
3141
3103
      ! contributions from others.
3142
3104
      ! NOTE: that the allocate interface assumes a contiguous