← 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-24 20:50:36 UTC
  • Revision ID: brendan.tollit05@imperial.ac.uk-20120724205036-gvzt2pnnm5c3h0gf
Add a default option for solver using GMRES with HYPRE preconditioner.

Adjust the dual porosity model to make absperm optional - turning this 
on now initiates the dual perm model. Adjust populate state as required 
and if no dual perm point the di variable to a constant field.

Show diffs side-by-side

added added

removed removed

Lines of Context:
main/Darcy_IMPES.F90
159
159
   type(darcy_impes_type) :: di
160
160
   type(darcy_impes_type) :: di_dual
161
161
   type(state_type) , dimension(:), pointer :: state_prime, state_dual
162
 
   logical :: have_dual, have_dual_pressure
 
162
   logical :: have_dual, have_dual_perm
163
163
   integer :: darcy_debug_log_unit, darcy_debug_err_unit
164
164
   integer :: number_phase_prime, number_phase_dual, number_of_first_adapts, p 
165
165
   character(len=OPTION_PATH_LEN) :: phase_name
256
256
   call set(time_field, current_time)
257
257
 
258
258
   ! ***** setting up dual *****
259
 
   have_dual =  have_option("/porous_media_dual")
260
 
   have_dual_pressure =  have_option("/porous_media_dual/have_dual_pressure")   
 
259
   have_dual      =  have_option("/porous_media_dual")
 
260
   have_dual_perm =  have_option("/porous_media_dual/scalar_field::AbsolutePermeabilityDual")   
261
261
   if (have_dual) then
262
262
      ! check there are the same number of prime and dual phases
263
263
      ! - NOTE: due to the options schemas all prime phases come first
287
287
                                  dt, &
288
288
                                  current_time, &
289
289
                                  have_dual, &
290
 
                                  have_dual_pressure, &
 
290
                                  have_dual_perm, &
291
291
                                  this_is_dual = .false.)
292
292
 
293
293
      call darcy_impes_initialise(di_dual, &
296
296
                                  dt, &
297
297
                                  current_time, &
298
298
                                  have_dual, &
299
 
                                  have_dual_pressure, &
 
299
                                  have_dual_perm, &
300
300
                                  this_is_dual = .true.)
301
301
   else
302
302
      ! *** Initialise data used in IMPES solver *** 
306
306
                                  dt, &
307
307
                                  current_time, &
308
308
                                  have_dual, &
309
 
                                  have_dual_pressure, &
 
309
                                  have_dual_perm, &
310
310
                                  this_is_dual = .false.)
311
311
   end if
312
312
   
507
507
   
508
508
   ! ***** Finalise dual permeability model *****
509
509
   if (have_dual) then
510
 
      call darcy_impes_finalise(di_dual, have_dual, have_dual_pressure, this_is_dual = .true.)
 
510
      call darcy_impes_finalise(di_dual, have_dual, have_dual_perm, this_is_dual = .true.)
511
511
   end if
512
512
   ! *** Finalise darcy impes variables ***
513
 
   call darcy_impes_finalise(di, have_dual, have_dual_pressure, this_is_dual = .false.)
 
513
   call darcy_impes_finalise(di, have_dual, have_dual_perm, this_is_dual = .false.)
514
514
    
515
515
   ! Deallocate state
516
516
   do i = 1, size(state)
564
564
                                     dt, &
565
565
                                     current_time, &
566
566
                                     have_dual, &
567
 
                                     have_dual_pressure, &
 
567
                                     have_dual_perm, &
568
568
                                     this_is_dual)
569
569
      
570
570
      !!< Initialise the Darcy IMPES type from options and state
575
575
      real,                                         intent(in)    :: dt
576
576
      real,                                         intent(in)    :: current_time
577
577
      logical ,                                     intent(in)    :: have_dual
578
 
      logical,                                      intent(in)    :: have_dual_pressure
 
578
      logical,                                      intent(in)    :: have_dual_perm
579
579
      logical ,                                     intent(in)    :: this_is_dual
580
580
      
581
581
      ! Local variables
605
605
      di%number_sele       = surface_element_count(di%pressure_mesh)
606
606
      di%number_pmesh_node = node_count(di%pressure_mesh)
607
607
      
 
608
      ! Allocate a field that is always zero on the element wise mesh to be pointed 
 
609
      ! at by capilliary pressure and saturation source if not required.
 
610
      ! NOTE using a FIELD_TYPE_CONSTANT causes issues.
 
611
      allocate(di%constant_zero_sfield_elementwisemesh)
 
612
      call allocate(di%constant_zero_sfield_elementwisemesh, di%elementwise_mesh)
 
613
      
 
614
      call zero(di%constant_zero_sfield_elementwisemesh)
 
615
      
608
616
      di%average_pressure         => extract_scalar_field(di%state(1), "AveragePressure")
609
 
      if (this_is_dual) then 
610
 
         di%porosity                     => extract_scalar_field(di%state(1), "PorosityDual")
611
 
         di%absolute_permeability        => extract_scalar_field(di%state(1), "AbsolutePermeabilityDual")
 
617
      if (have_dual .and. this_is_dual) then 
 
618
         di%porosity              => extract_scalar_field(di%state(1), "PorosityDual")
 
619
         if (have_dual_perm) then
 
620
            di%absolute_permeability => extract_scalar_field(di%state(1), "AbsolutePermeabilityDual")
 
621
         else
 
622
            di%absolute_permeability => di%constant_zero_sfield_elementwisemesh
 
623
         end if
612
624
      else
613
625
         di%porosity              => extract_scalar_field(di%state(1), "Porosity")
614
626
         di%absolute_permeability => extract_scalar_field(di%state(1), "AbsolutePermeability")
680
692
      call allocate(di%matrix, di%sparsity_pmesh_pmesh)
681
693
      ! Only allocate the pressure matrix for the prime di call here
682
694
      if (.not. this_is_dual) then
683
 
         if (have_dual .and. have_dual_pressure) then
 
695
         if (have_dual .and. have_dual_perm) then
684
696
            call allocate(di%dual_block_pressure_matrix, di%sparsity_pmesh_pmesh, blocks=(/2,2/), name ='DualBlockPressureMatrix')
685
697
         else
686
698
            call allocate(di%pressure_matrix, di%sparsity_pmesh_pmesh, name ='PressureMatrix')      
1432
1444
 
1433
1445
! ----------------------------------------------------------------------------
1434
1446
   
1435
 
   subroutine darcy_impes_finalise(di, have_dual, have_dual_pressure, this_is_dual)
 
1447
   subroutine darcy_impes_finalise(di, have_dual, have_dual_perm, this_is_dual)
1436
1448
      
1437
1449
      !!< Finalise (ie deallocate) the Darcy IMPES data
1438
1450
      
1439
1451
      type(darcy_impes_type), intent(inout) :: di
1440
1452
      logical,                intent(in)    :: have_dual
1441
 
      logical,                intent(in)    :: have_dual_pressure
 
1453
      logical,                intent(in)    :: have_dual_perm
1442
1454
      logical,                intent(in)    :: this_is_dual
1443
1455
          
1444
1456
      ! Deallocate, nullify or zero Darcy IMPES data
1473
1485
      di%number_sele       = 0
1474
1486
      di%number_pmesh_node = 0
1475
1487
      
 
1488
      call deallocate(di%constant_zero_sfield_elementwisemesh)
 
1489
      deallocate(di%constant_zero_sfield_elementwisemesh)
 
1490
      
1476
1491
      nullify(di%average_pressure)      
1477
1492
      nullify(di%porosity)
1478
1493
      nullify(di%absolute_permeability)
1504
1519
      nullify(di%sparsity_pmesh_pmesh)
1505
1520
      call deallocate(di%matrix)
1506
1521
      if (.not. this_is_dual) then 
1507
 
         if (have_dual .and. have_dual_pressure) then
 
1522
         if (have_dual .and. have_dual_perm) then
1508
1523
            call deallocate(di%dual_block_pressure_matrix)
1509
1524
         else
1510
1525
            call deallocate(di%pressure_matrix)     
1697
1712
                                                    dt, &
1698
1713
                                                    current_time, &
1699
1714
                                                    have_dual, &
1700
 
                                                    have_dual_pressure, &
 
1715
                                                    have_dual_perm, &
1701
1716
                                                    this_is_dual)
1702
1717
      
1703
1718
      !!< Update the Darcy IMPES data post spatial adapt
1708
1723
      real,                                         intent(in)    :: dt
1709
1724
      real,                                         intent(in)    :: current_time
1710
1725
      logical,                                      intent(in)    :: have_dual  
1711
 
      logical,                                      intent(in)    :: have_dual_pressure
 
1726
      logical,                                      intent(in)    :: have_dual_perm
1712
1727
      logical,                                      intent(in)    :: this_is_dual  
1713
1728
        
1714
1729
      ewrite(1,*) 'Update Darcy IMPES data post spatial adapt'
1716
1731
      ! ALL THE IMPORTANT SOLUTION DATA IS IN STATE
1717
1732
      ! WHICH IS NOT DEALLOCATED IN THE FOLLOWING PROCEDURE
1718
1733
      
1719
 
      call darcy_impes_finalise(di, have_dual, have_dual_pressure, this_is_dual)
 
1734
      call darcy_impes_finalise(di, have_dual, have_dual_perm, this_is_dual)
1720
1735
            
1721
1736
      call darcy_impes_initialise(di, &
1722
1737
                                  state, &
1724
1739
                                  dt, &
1725
1740
                                  current_time, &
1726
1741
                                  have_dual, &
1727
 
                                  have_dual_pressure, &
 
1742
                                  have_dual_perm, &
1728
1743
                                  this_is_dual)
1729
1744
      
1730
1745
      ewrite(1,*) 'Finished updating Darcy IMPES data post spatial adapt'
1770
1785
                                                    dt, &
1771
1786
                                                    current_time, &
1772
1787
                                                    have_dual, &
1773
 
                                                    have_dual_pressure, &
 
1788
                                                    have_dual_perm, &
1774
1789
                                                    this_is_dual = .false.)
1775
1790
 
1776
1791
         call darcy_impes_update_post_spatial_adapt(di_dual, &
1779
1794
                                                    dt, &
1780
1795
                                                    current_time, &
1781
1796
                                                    have_dual, &
1782
 
                                                    have_dual_pressure, &
 
1797
                                                    have_dual_perm, &
1783
1798
                                                    this_is_dual = .true.)
1784
1799
      else                                                                 
1785
1800
         ! *** Update Darcy IMPES post spatial adapt ***
1789
1804
                                                    dt, &
1790
1805
                                                    current_time, &
1791
1806
                                                    have_dual, &
1792
 
                                                    have_dual_pressure, &
 
1807
                                                    have_dual_perm, &
1793
1808
                                                    this_is_dual = .false.)
1794
1809
      end if  
1795
1810
 
1907
1922
         if (have_dual) call darcy_impes_assemble_and_solve_part_one(di_dual, have_dual)
1908
1923
         
1909
1924
         ! This one solves for the pressure phase 1
1910
 
         call darcy_impes_assemble_and_solve_part_two(di, di_dual, have_dual, have_dual_pressure)
 
1925
         call darcy_impes_assemble_and_solve_part_two(di, di_dual, have_dual, have_dual_perm)
1911
1926
 
1912
1927
         call darcy_impes_assemble_and_solve_part_three(di, have_dual)
1913
1928
         if (have_dual) call darcy_impes_assemble_and_solve_part_three(di_dual, have_dual)