Commits

Guoping Tang committed 8ac7fc2

Replace the regulating function with Nathan's f(x) = (1 - x^2)^2 or 1 - (1 - x^2)^2 as Gautam suggested for source/sink downregulation.

Comments (0)

Files changed (3)

regression_tests/ngee/sink.regression.gold

 -- PRESSURE: Liquid Pressure --
-      Max:  -5.5626002228507E+04
-      Min:  -4.5479426188635E+05
-     Mean:  -1.0515629012847E+05
-        1:  -5.5626002228507E+04
-       10:  -6.4512724399693E+04
-       19:  -7.4860774716110E+04
-       28:  -9.9552326363003E+04
+      Max:  -5.5626963027381E+04
+      Min:  -1.0996011388029E+06
+     Mean:  -1.3412263717370E+05
+        1:  -5.5626963027381E+04
+       10:  -6.4530359665847E+04
+       19:  -7.5178264233936E+04
+       28:  -1.0346647838158E+05
 -- SATURATION: Liquid Saturation --
-      Max:   3.6084799701071E-01
-      Min:   2.5019971345394E-01
-     Mean:   3.3855219658553E-01
-        1:   3.6084799701071E-01
-       10:   3.5498831480167E-01
-       19:   3.4866969719247E-01
-       28:   3.3540519502187E-01
+      Max:   3.6084733981985E-01
+      Min:   2.0296649715393E-01
+     Mean:   3.3462616044583E-01
+        1:   3.6084733981985E-01
+       10:   3.5497710373391E-01
+       19:   3.4848368064111E-01
+       28:   3.3350188618848E-01
 -- DISCRETE: Material ID --
       Max:         1
       Min:         1
        19:         1
        28:         1
 -- SOLUTION: Flow --
-   Time (seconds):   1.7262787818909E+00
+   Time (seconds):   2.2709729671478E+00
    Time Steps:          381
-   Newton Iterations:          816
-   Solver Iterations:          816
+   Newton Iterations:         1118
+   Solver Iterations:         1118
    Time Step Cuts:            0
-   Solution 2-Norm:   8.0882976468920E+05
-   Residual 2-Norm:   1.4772022745979E-18
+   Solution 2-Norm:   1.3914910224507E+06
+   Residual 2-Norm:   1.3619141086517E-18

regression_tests/ngee/source.regression.gold

 -- PRESSURE: Liquid Pressure --
-      Max:   6.7273529652321E+06
-      Min:   6.6910643622375E+06
-     Mean:   6.7092085732808E+06
-        1:   6.7273529652321E+06
-       10:   6.7185259050700E+06
-       19:   6.7096989109137E+06
-       28:   6.7008719827643E+06
+      Max:   9.3986952793657E+05
+      Min:   9.0367406812395E+05
+     Mean:   9.2177175035915E+05
+        1:   9.3986952793657E+05
+       10:   9.3106517278821E+05
+       19:   9.2226085242631E+05
+       28:   9.1345656685140E+05
 -- SATURATION: Liquid Saturation --
       Max:   1.0000000000000E+00
       Min:   1.0000000000000E+00
        19:         1
        28:         1
 -- SOLUTION: Flow --
-   Time (seconds):   2.1375660896301E+00
-   Time Steps:          381
-   Newton Iterations:         1126
-   Solver Iterations:         1126
-   Time Step Cuts:            0
-   Solution 2-Norm:   4.1358392416848E+07
-   Residual 2-Norm:   1.8034576377192E-15
+   Time (seconds):   1.7880368232727E+00
+   Time Steps:          389
+   Newton Iterations:          906
+   Solver Iterations:          906
+   Time Step Cuts:            3
+   Solution 2-Norm:   5.6825674663828E+06
+   Residual 2-Norm:   2.7164624746046E-16

src/pflotran/srcsink_sandbox_downreg.F90

 ! Source/sink Sandbox for downregulating source or sink terms to avoid 
 ! overpressurization (too big a pressure, + or -)
 ! the source is regulated by multiplying it by 
-!   pressure_max/(pressure + pressure_max)
+!   (1 - x^2)^2 
+! with 
+!   x = (pressure - pref - pressure_max - delta/2)/delta
+! for pressure in between of x*delta and (x+1)*delta
 ! the sink is regulated by multiplying it by 
-!   pressure_min/(pressure + pressure_min)
+!   1 - (1 - x^2)^2 
+! with 
+!   x = (pressure - pref - pressure_min - delta/2)/delta
+! for pressure in between of x*delta and (x+1)*delta
+! for transpiration sink, the physical interpretation can be that the soil water
+! extration rate is decreased from q to 0 as pressure decreases from 
+! pressure_min - delta/2 to pressure_min + delta/2
+
 ! extended from srcsink_sandbox_mass_rate
 ! currently work in RICHARDS mode
   
     PetscReal, pointer :: rate(:)
     PetscReal :: pressure_min       ! Pa
     PetscReal :: pressure_max       ! Pa
+    PetscReal :: pressure_delta     ! Pa
   contains
     procedure, public :: ReadInput => DownregRead
     procedure, public :: Setup => DownregSetup
   nullify(DownregCreate%rate)
   DownregCreate%pressure_min = -1.0d9
   DownregCreate%pressure_max = 1.0d9
+  DownregCreate%pressure_delta = 1.0d4
   
 end function DownregCreate
 
         call InputReadDouble(input,option,this%pressure_max)
         call InputErrorMsg(input,option,'maximum pressure (Pa)', &
           'SOURCE_SINK_SANDBOX,DOWNREG')
-
       case('NEGATIVE_REG_PRESSURE')
         call InputReadDouble(input,option,this%pressure_min)
         call InputErrorMsg(input,option,'minimum pressure (Pa)', &
           'SOURCE_SINK_SANDBOX,DOWNREG')
+      case('DELTA_REG_PRESSURE')
+        call InputReadDouble(input,option,this%pressure_delta)
+        call InputErrorMsg(input,option,'delta pressure (Pa)', &
+          'SOURCE_SINK_SANDBOX,DOWNREG')
+        if (this%pressure_delta <= 1.0d-10) then
+          option%io_buffer = 'SRCSINK_SANDBOX,DOWNREG,DELTA_REG_PRESSURE' // &
+            ': the pressure delta is too close to 0 Pa.'
+          call printErrMsg(option)
+        endif 
       case default
         option%io_buffer = 'SRCSINK_SANDBOX,DOWNREG keyword: ' // &
           trim(word) // ' not recognized.'
   ! Author: Guoping Tang
   ! Date: 06/04/14
   ! Glenn suggested to use reference pressure for RICHARDS mode 
+  ! Gautam suggested to use Nathan's smooth function (6/8/2014)
 
   use Option_module
   use Reaction_Aux_module
   class(material_auxvar_type) :: material_auxvar
   PetscReal :: aux_real(:)
   PetscReal :: pressure
+  PetscReal :: pressure_lower, pressure_upper, x
   PetscReal :: rate_regulator  
   PetscReal :: drate_regulator  
   PetscReal :: temp_real
       ! regulate liquid pressure in Richards mode
       pressure = aux_real(1)
       if (this%rate(idof) > 0.0d0) then
-        ! source (try to avoid this%pressure_max + pressure = 0)
-        ! regulate the source under wet (not dry) condition 
-        ! not the case with too much infiltration to dry soil
-        if (pressure > option%reference_pressure) then
-          rate_regulator = this%pressure_max / &
-            (this%pressure_max + pressure - option%reference_pressure)
+        ! source
+        pressure_lower = this%pressure_max - this%pressure_delta/2.0d0 &
+                                           - option%reference_pressure
+        pressure_upper = pressure_lower + this%pressure_delta
+        if (pressure <= pressure_lower) then
+          rate_regulator = 1.0d0
+          drate_regulator = 0.0d0
+        elseif (pressure >= pressure_upper) then
+          rate_regulator = 0.0d0
+          drate_regulator = 0.0d0
         else
-          rate_regulator = 1.0d0 
+          x = (pressure - pressure_lower)/this%pressure_delta
+          rate_regulator = (1.0d0 - x * x) * (1.0d0 - x * x) 
+          drate_regulator = 2.0d0 * (1.0d0 - x * x) * (-2.0d0) * x / &
+            this%pressure_delta
         endif
-          Residual(idof) = this%rate(idof) * rate_regulator
+        Residual(idof) = this%rate(idof) * rate_regulator
       else
-        ! sink (try to avoid this%pressure_min + pressure = 0)
-        ! regulate the sink under dry (not wet) condition, 
-        if (pressure < option%reference_pressure) then
-          rate_regulator = this%pressure_min / &
-            (this%pressure_min + pressure - option%reference_pressure)
+        ! sink
+        pressure_lower = this%pressure_min - this%pressure_delta/2.0d0 &
+                                           - option%reference_pressure
+        pressure_upper = pressure_lower + this%pressure_delta
+        if (pressure <= pressure_lower) then
+          rate_regulator = 0.0d0
+          drate_regulator = 0.0d0
+        elseif (pressure >= pressure_upper) then
+          rate_regulator = 1.0d0
+          drate_regulator = 0.0d0
         else
-          rate_regulator = 1.0d0
+          x = (pressure - pressure_lower)/this%pressure_delta
+          rate_regulator = 1.0d0 - (1.0d0 - x * x) * (1.0d0 - x * x) 
+          drate_regulator = (-2.0d0) * (1.0d0 - x * x) * (-2.0d0) * x / &
+            this%pressure_delta
         endif
         Residual(idof) = this%rate(idof) * rate_regulator
       endif
     do idof = 1, option%nflowdof
       if (option%iflowmode == RICHARDS_MODE .and. idof == ONE_INTEGER) then
         ! regulate liquid pressure in Richards mode
-        pressure = aux_real(1)
-        if (this%rate(idof) > 0.0d0) then
-          ! source
-          if (pressure > option%reference_pressure) then
-            temp_real = this%pressure_max - option%reference_pressure 
-            drate_regulator = (-1.0d0) * this%pressure_max / &
-              (temp_real + pressure) / (temp_real + pressure)
-          else
-            drate_regulator = 0.0d0
-          endif
-
-          Jacobian(idof,idof) = -1.0d0 * this%rate(idof) * drate_regulator
-        else
-          ! sink           
-          if (pressure < option%reference_pressure) then
-            temp_real = this%pressure_min - option%reference_pressure 
-            drate_regulator = -1.0d0 * this%pressure_min / &
-              (temp_real + pressure) / (temp_real + pressure)
-          else
-            drate_regulator = 0.0d0
-          endif
-          Jacobian(idof,idof) = -1.0d0 * this%rate(idof) * drate_regulator
-        endif
+        Jacobian(idof,idof) = -1.0d0 * this%rate(idof) * drate_regulator
       else
         ! since the rates are constant, there is no derivative
       endif