← Back to branch summary

~fluidity-core/fluidity/darcy_weak_bcs

« back to all changes in this revision

Viewing changes to main/Darcy_IMPES.F90

  • Committer: Brendan Tollit
  • Date: 2012-07-19 13:58:00 UTC
  • Revision ID: brendan.tollit05@imperial.ac.uk-20120719135800-6u8xid86jmgu0222
Add an option have_dual_pressure.

Correct a couple of bugs for the dual model.

Only use the block pressure matrix when really required as it is 
bugged in 2d.

Show diffs side-by-side

added added

removed removed

Lines of Context:
main/Darcy_IMPES.F90
156
156
   type(darcy_impes_type) :: di
157
157
   type(darcy_impes_type) :: di_dual
158
158
   type(state_type) , dimension(:), pointer :: state_prime, state_dual
159
 
   logical :: have_dual
 
159
   logical :: have_dual, have_dual_pressure
160
160
   integer :: darcy_debug_log_unit, darcy_debug_err_unit   
161
161
   
162
162
#ifdef HAVE_ZOLTAN
238
238
 
239
239
   ! ***** setting up dual *****
240
240
   have_dual =  have_option("/porous_media_dual")
 
241
   have_dual_pressure =  have_option("/porous_media_dual/have_dual_pressure")   
241
242
   if (have_dual) then
242
243
      ! this is set up for two phases in each porous media
243
244
      if ( option_count("/material_phase") /= 4 ) then
250
251
                                  dt, &
251
252
                                  current_time, &
252
253
                                  have_dual, &
 
254
                                  have_dual_pressure, &
253
255
                                  this_is_dual = .false.)
254
256
 
255
257
      call darcy_impes_initialise(di_dual, &
257
259
                                  dt, &
258
260
                                  current_time, &
259
261
                                  have_dual, &
 
262
                                  have_dual_pressure, &
260
263
                                  this_is_dual = .true.)
261
264
   else
262
265
      ! *** Initialise data used in IMPES solver *** 
265
268
                                  dt, &
266
269
                                  current_time, &
267
270
                                  have_dual, &
 
271
                                  have_dual_pressure, &
268
272
                                  this_is_dual = .false.)
269
273
   end if
270
274
   
427
431
         call darcy_impes_assemble_and_solve_part_one(di, have_dual)
428
432
         if (have_dual) call darcy_impes_assemble_and_solve_part_one(di_dual, have_dual)
429
433
         
430
 
         call darcy_impes_assemble_and_solve_part_two(di, di_dual, have_dual)
 
434
         ! This one solves for the pressure phase 1
 
435
         call darcy_impes_assemble_and_solve_part_two(di, di_dual, have_dual, have_dual_pressure)
431
436
 
432
437
         call darcy_impes_assemble_and_solve_part_three(di, have_dual)
433
438
         if (have_dual) call darcy_impes_assemble_and_solve_part_three(di_dual, have_dual)
545
550
                                                       dt, &
546
551
                                                       current_time, &
547
552
                                                       have_dual, &
 
553
                                                       have_dual_pressure, &
548
554
                                                       this_is_dual = .false.)
549
555
                                                       
550
556
               call darcy_impes_update_post_spatial_adapt(di_dual, &
552
558
                                                       dt, &
553
559
                                                       current_time, &
554
560
                                                       have_dual, &
 
561
                                                       have_dual_pressure, &
555
562
                                                       this_is_dual = .true.)
556
563
            else                                                                 
557
564
               ! *** Update Darcy IMPES post spatial adapt ***
560
567
                                                       dt, &
561
568
                                                       current_time, &
562
569
                                                       have_dual, &
 
570
                                                       have_dual_pressure, &
563
571
                                                       this_is_dual = .false.)
564
572
            end if  
565
573
                                                                            
635
643
   
636
644
   ! ***** Finalise dual permeability model *****
637
645
   if (have_dual) then
638
 
      call darcy_impes_finalise(di_dual, have_dual)
 
646
      call darcy_impes_finalise(di_dual, have_dual, have_dual_pressure, this_is_dual = .true.)
639
647
   end if
640
648
   ! *** Finalise darcy impes variables ***
641
 
   call darcy_impes_finalise(di, have_dual)
 
649
   call darcy_impes_finalise(di, have_dual, have_dual_pressure, this_is_dual = .false.)
642
650
    
643
651
   ! Deallocate state
644
652
   do i = 1, size(state)
691
699
                                     dt, &
692
700
                                     current_time, &
693
701
                                     have_dual, &
 
702
                                     have_dual_pressure, &
694
703
                                     this_is_dual)
695
704
      
696
705
      !!< Initialise the Darcy IMPES type from options and state
700
709
      real,                                         intent(in)    :: dt
701
710
      real,                                         intent(in)    :: current_time
702
711
      logical ,                                     intent(in)    :: have_dual
 
712
      logical,                                      intent(in)    :: have_dual_pressure
703
713
      logical ,                                     intent(in)    :: this_is_dual
704
714
      
705
715
      ! Local variables
810
820
      call allocate(di%matrix, di%sparsity_pmesh_pmesh)
811
821
      ! Only allocate the pressure matrix for the prime di call here
812
822
      if (.not. this_is_dual) then
813
 
         if (have_dual) then
814
 
            call allocate(di%pressure_matrix, di%sparsity_pmesh_pmesh, blocks=(/2,2/), name ='PressureMatrix')
 
823
         if (have_dual .and. have_dual_pressure) then
 
824
            call allocate(di%dual_block_pressure_matrix, di%sparsity_pmesh_pmesh, blocks=(/2,2/), name ='DualBlockPressureMatrix')
815
825
         else
816
 
            call allocate(di%pressure_matrix, di%sparsity_pmesh_pmesh, blocks=(/1,1/), name ='PressureMatrix')      
 
826
            call allocate(di%pressure_matrix, di%sparsity_pmesh_pmesh, name ='PressureMatrix')      
817
827
         end if
818
828
      end if
819
829
      allocate(di%pressure_rhs)
1554
1564
 
1555
1565
! ----------------------------------------------------------------------------
1556
1566
   
1557
 
   subroutine darcy_impes_finalise(di, have_dual)
 
1567
   subroutine darcy_impes_finalise(di, have_dual, have_dual_pressure, this_is_dual)
1558
1568
      
1559
1569
      !!< Finalise (ie deallocate) the Darcy IMPES data
1560
1570
      
1561
1571
      type(darcy_impes_type), intent(inout) :: di
1562
1572
      logical,                intent(in)    :: have_dual
1563
 
            
 
1573
      logical,                intent(in)    :: have_dual_pressure
 
1574
      logical,                intent(in)    :: this_is_dual
 
1575
          
1564
1576
      ! Deallocate, nullify or zero Darcy IMPES data
1565
1577
      !  - pointers to data in state are nullified
1566
1578
      !  - objects with memory only in darcy_impes_type are deallocated
1625
1637
      call deallocate(di%positions_pressure_mesh)
1626
1638
      nullify(di%sparsity_pmesh_pmesh)
1627
1639
      call deallocate(di%matrix)
1628
 
      call deallocate(di%pressure_matrix)
 
1640
      if (.not. this_is_dual) then 
 
1641
         if (have_dual .and. have_dual_pressure) then
 
1642
            call deallocate(di%dual_block_pressure_matrix)
 
1643
         else
 
1644
            call deallocate(di%pressure_matrix)     
 
1645
         end if
 
1646
      end if
1629
1647
      call deallocate(di%pressure_rhs)
1630
1648
      deallocate(di%pressure_rhs)
1631
1649
      call deallocate(di%lhs)
1808
1826
                                                    dt, &
1809
1827
                                                    current_time, &
1810
1828
                                                    have_dual, &
 
1829
                                                    have_dual_pressure, &
1811
1830
                                                    this_is_dual)
1812
1831
      
1813
1832
      !!< Update the Darcy IMPES data post spatial adapt
1817
1836
      real,                                         intent(in)    :: dt
1818
1837
      real,                                         intent(in)    :: current_time
1819
1838
      logical,                                      intent(in)    :: have_dual  
 
1839
      logical,                                      intent(in)    :: have_dual_pressure
1820
1840
      logical,                                      intent(in)    :: this_is_dual  
1821
1841
        
1822
1842
      ewrite(1,*) 'Update Darcy IMPES data post spatial adapt'
1824
1844
      ! ALL THE IMPORTANT SOLUTION DATA IS IN STATE
1825
1845
      ! WHICH IS NOT DEALLOCATED IN THE FOLLOWING PROCEDURE
1826
1846
      
1827
 
      call darcy_impes_finalise(di, have_dual)
 
1847
      call darcy_impes_finalise(di, have_dual, have_dual_pressure, this_is_dual)
1828
1848
            
1829
1849
      call darcy_impes_initialise(di, &
1830
1850
                                  state, &
1831
1851
                                  dt, &
1832
1852
                                  current_time, &
1833
 
                                  this_is_dual, &
1834
 
                                  have_dual)
 
1853
                                  have_dual, &
 
1854
                                  have_dual_pressure, &
 
1855
                                  this_is_dual)
1835
1856
      
1836
1857
      ewrite(1,*) 'Finished updating Darcy IMPES data post spatial adapt'
1837
1858