Commits

Anonymous committed 239ac55

Re-arranged option%flow_dt in residual calculations in TH, THC, THMC.

Comments (0)

Files changed (4)

regression_tests/default/anisothermal/thc_1d.regression.gold

 -- GENERIC: Temperature --
-      Max:   9.6760222587887E+01
-      Min:   2.5001659582776E+01
-     Mean:   3.3727078188315E+01
-        1:   2.5001659582776E+01
-       51:   2.6128270179649E+01
-      100:   9.6760222587887E+01
-        1:   2.5001659582776E+01
-       51:   2.6128270179649E+01
+      Max:   9.6760222627490E+01
+      Min:   2.5001659583837E+01
+     Mean:   3.3727078619022E+01
+        1:   2.5001659583837E+01
+       51:   2.6128270699353E+01
+      100:   9.6760222627490E+01
+        1:   2.5001659583837E+01
+       51:   2.6128270699353E+01
 -- PRESSURE: Liquid Pressure --
-      Max:   9.9517919076094E+05
-      Min:   1.5162788207926E+05
-     Mean:   5.3458322001463E+05
-        1:   9.9517919076094E+05
-       51:   5.1590008309070E+05
-      100:   1.5162788207926E+05
-        1:   9.9517919076094E+05
-       51:   5.1590008309070E+05
+      Max:   9.9517919071995E+05
+      Min:   1.5162788209189E+05
+     Mean:   5.3458321881632E+05
+        1:   9.9517919071995E+05
+       51:   5.1590008084006E+05
+      100:   1.5162788209189E+05
+        1:   9.9517919071995E+05
+       51:   5.1590008084006E+05
 -- SATURATION: Liquid Saturation --
       Max:   1.0000000000000E+00
       Min:   1.0000000000000E+00
         1:   1.0000000000000E+00
        51:   1.0000000000000E+00
 -- GENERIC: Liquid Density --
-      Max:   9.9756033956781E+02
-      Min:   9.6046540737091E+02
-     Mean:   9.9385348480257E+02
-        1:   9.9756033956781E+02
-       51:   9.9704854156178E+02
-      100:   9.6046540737091E+02
-        1:   9.9756033956781E+02
-       51:   9.9704854156178E+02
+      Max:   9.9756033956752E+02
+      Min:   9.6046540734286E+02
+     Mean:   9.9385348465539E+02
+        1:   9.9756033956752E+02
+       51:   9.9704854142089E+02
+      100:   9.6046540734286E+02
+        1:   9.9756033956752E+02
+       51:   9.9704854142089E+02
 -- GENERIC: Liquid Energy --
-      Max:   7.3014959477002E+00
-      Min:   1.8860968863196E+00
-     Mean:   2.5439927286570E+00
-        1:   1.8860968863196E+00
-       51:   1.9715607346130E+00
-      100:   7.3014959477002E+00
-        1:   1.8860968863196E+00
-       51:   1.9715607346130E+00
+      Max:   7.3014959507125E+00
+      Min:   1.8860968863961E+00
+     Mean:   2.5439927610851E+00
+        1:   1.8860968863961E+00
+       51:   1.9715607737319E+00
+      100:   7.3014959507125E+00
+        1:   1.8860968863961E+00
+       51:   1.9715607737319E+00
 -- GENERIC: Liquid Viscosity --
-      Max:   8.9034195495354E-04
-      Min:   2.8878020354658E-04
-     Mean:   7.8011545875108E-04
-        1:   8.9034195495354E-04
-       51:   8.6794400989330E-04
-      100:   2.8878020354658E-04
-        1:   8.9034195495354E-04
-       51:   8.6794400989330E-04
+      Max:   8.9034195493200E-04
+      Min:   2.8878020342311E-04
+     Mean:   7.8011545204717E-04
+        1:   8.9034195493200E-04
+       51:   8.6794399975086E-04
+      100:   2.8878020342311E-04
+        1:   8.9034195493200E-04
+       51:   8.6794399975086E-04
 -- GENERIC: Liquid Mobility --
-      Max:   3.4628412464523E+03
-      Min:   1.1231639646277E+03
-     Mean:   1.3815239281469E+03
-        1:   1.1231639646277E+03
-       51:   1.1521480517193E+03
-      100:   3.4628412464523E+03
-        1:   1.1231639646277E+03
-       51:   1.1521480517193E+03
+      Max:   3.4628412479329E+03
+      Min:   1.1231639646549E+03
+     Mean:   1.3815239401681E+03
+        1:   1.1231639646549E+03
+       51:   1.1521480651828E+03
+      100:   3.4628412479329E+03
+        1:   1.1231639646549E+03
+       51:   1.1521480651828E+03
 -- GENERIC: Liquid Mole Fraction H2O --
       Max:   1.0000000000000E+00
       Min:   1.0000000000000E+00
         1:   0.0000000000000E+00
        51:   0.0000000000000E+00
 -- GENERIC: pH --
-      Max:   7.3589063696322E+00
-      Min:   5.0421807971060E+00
-     Mean:   7.2357962646775E+00
-        1:   5.0421807971060E+00
-       51:   7.3441569729744E+00
-      100:   6.6652550497889E+00
-        1:   5.0421807971060E+00
-       51:   7.3441569729744E+00
+      Max:   7.3589063695209E+00
+      Min:   5.0421807965375E+00
+     Mean:   7.2357962587363E+00
+        1:   5.0421807965375E+00
+       51:   7.3441569657174E+00
+      100:   6.6652550481804E+00
+        1:   5.0421807965375E+00
+       51:   7.3441569657174E+00
 -- CONCENTRATION: Total H+ --
-      Max:   1.4773810851462E-03
-      Min:   2.1524203235240E-04
-     Mean:   2.7330855582393E-04
-        1:   1.4773810851462E-03
-       51:   2.1891676747933E-04
-      100:   6.9599216549039E-04
-        1:   1.4773810851462E-03
-       51:   2.1891676747933E-04
+      Max:   1.4773810852009E-03
+      Min:   2.1524203238062E-04
+     Mean:   2.7330855849252E-04
+        1:   1.4773810852009E-03
+       51:   2.1891676959690E-04
+      100:   6.9599216951757E-04
+        1:   1.4773810852009E-03
+       51:   2.1891676959690E-04
 -- CONCENTRATION: Total HCO3- --
-      Max:   2.8070386416937E-03
-      Min:   1.5449485058711E-03
-     Mean:   2.7333994735882E-03
-        1:   1.5449485058711E-03
-       51:   2.8014606928362E-03
-      100:   2.2069807690324E-03
-        1:   1.5449485058711E-03
-       51:   2.8014606928362E-03
+      Max:   2.8070386416645E-03
+      Min:   1.5449485058160E-03
+     Mean:   2.7333994740141E-03
+        1:   1.5449485058160E-03
+       51:   2.8014606923692E-03
+      100:   2.2069807754805E-03
+        1:   1.5449485058160E-03
+       51:   2.8014606923692E-03
 -- CONCENTRATION: Total Ca++ --
-      Max:   1.8080578423991E-03
-      Min:   5.4595153777427E-04
-     Mean:   1.7395659639201E-03
-        1:   5.4595153777427E-04
-       51:   1.8031089363802E-03
-      100:   1.2474355799222E-03
-        1:   5.4595153777427E-04
-       51:   1.8031089363802E-03
+      Max:   1.8080578423688E-03
+      Min:   5.4595153771912E-04
+     Mean:   1.7395659627193E-03
+        1:   5.4595153771912E-04
+       51:   1.8031089350086E-03
+      100:   1.2474355810951E-03
+        1:   5.4595153771912E-04
+       51:   1.8031089350086E-03
 -- VOLUME_FRACTION: Calcite VF --
-      Max:   1.0259970533119E-01
+      Max:   1.0260077226923E-01
       Min:   0.0000000000000E+00
-     Mean:   9.8740095856125E-02
+     Mean:   9.8740260620933E-02
         1:   0.0000000000000E+00
-       51:   1.0000699802974E-01
-      100:   1.0259970533119E-01
+       51:   1.0000700260352E-01
+      100:   1.0260077226923E-01
         1:   0.0000000000000E+00
-       51:   1.0000699802974E-01
+       51:   1.0000700260352E-01
 -- RATE: Calcite Rate --
-      Max:   2.5097340432942E-09
-      Min:  -5.6618101532350E-08
-     Mean:  -5.1544899493714E-10
+      Max:   2.5097335585217E-09
+      Min:  -5.6618102004787E-08
+     Mean:  -5.1544915160033E-10
         1:   0.0000000000000E+00
-       51:   2.3361179718785E-12
-      100:   2.5097340432942E-09
+       51:   2.3360947344598E-12
+      100:   2.5097335585217E-09
         1:   0.0000000000000E+00
-       51:   2.3361179718785E-12
+       51:   2.3360947344598E-12
 -- DISCRETE: Material ID --
       Max:         1
       Min:         1
         1:         1
        51:         1
 -- GENERIC: LIQUID VELOCITY [m/y] --
-        1:   3.4149516081499E-01  0.0000000000000E+00  0.0000000000000E+00
-       51:   3.4167063584128E-01  0.0000000000000E+00  0.0000000000000E+00
-      100:   3.5473330958586E-01  0.0000000000000E+00  0.0000000000000E+00
-        1:   3.4149516081499E-01  0.0000000000000E+00  0.0000000000000E+00
-       51:   3.4167063584128E-01  0.0000000000000E+00  0.0000000000000E+00
+        1:   3.4149516371876E-01  0.0000000000000E+00  0.0000000000000E+00
+       51:   3.4167063877721E-01  0.0000000000000E+00  0.0000000000000E+00
+      100:   3.5473331250039E-01  0.0000000000000E+00  0.0000000000000E+00
+        1:   3.4149516371876E-01  0.0000000000000E+00  0.0000000000000E+00
+       51:   3.4167063877721E-01  0.0000000000000E+00  0.0000000000000E+00
 -- SOLUTION: Flow --
-   Time (seconds):   2.0868527889252E+00
+   Time (seconds):   9.7909355163574E-01
    Time Steps:          172
-   Newton Iterations:          390
-   Solver Iterations:          390
+   Newton Iterations:          211
+   Solver Iterations:          211
    Time Step Cuts:            0
-   Solution 2-Norm:   5.9342757408132E+06
-   Residual 2-Norm:   2.5324670047437E-09
+   Solution 2-Norm:   5.9342757286392E+06
+   Residual 2-Norm:   4.0331334175962E-17
 -- SOLUTION: Transport --
-   Time (seconds):   4.6810803413391E+00
+   Time (seconds):   5.0090203285217E+00
    Time Steps:          173
    Newton Iterations:          656
    Solver Iterations:          656
    Time Step Cuts:            0
-   Solution 2-Norm:   2.9787114317876E-02
-   Residual 2-Norm:   5.6630814576365E-14
+   Solution 2-Norm:   2.9787114296342E-02
+   Residual 2-Norm:   5.0905228878917E-14

src/pflotran/th.F90

 #endif                    
 #endif
 
-
+  J = J/option%flow_dt
   J(option%nflowdof,:) = vol_frac_prim*J(option%nflowdof,:)
 
   if (option%numerical_derivatives_flow) then
 #endif
 #endif
 
-  Res(1:option%nflowdof-1) = mol(:)
-  Res(option%nflowdof) = vol_frac_prim*eng
+  Res(1:option%nflowdof-1) = mol(:)/option%flow_dt
+  Res(option%nflowdof) = vol_frac_prim*eng/option%flow_dt
 
 end subroutine THAccumulation
 
   Jdn(option%nflowdof,2) = Jdn(option%nflowdof,2) + Dk*area*(-1.d0) + &
                            area*(global_aux_var_up%temp(1) - & 
                            global_aux_var_dn%temp(1))*dDk_dt_dn 
-  Jup = Jup*option%flow_dt
-  Jdn = Jdn*option%flow_dt
+
  ! note: Res is the flux contribution, for node up J = J + Jup
  !                                              dn J = J - Jdn  
 
   cond = Dk*area*(global_aux_var_up%temp(1) - global_aux_var_dn%temp(1)) 
   fluxe = fluxe + cond
 
-  Res(1:option%nflowdof-1) = fluxm(:) * option%flow_dt
-  Res(option%nflowdof) = fluxe * option%flow_dt
+  Res(1:option%nflowdof-1) = fluxm(:)
+  Res(option%nflowdof) = fluxe
   
  ! note: Res is the flux contribution, for node 1 R = R + Res_FL
  !                                              2 R = R - Res_FL  
 
   end select
 
-  Jdn = Jdn * option%flow_dt
-
   if (option%numerical_derivatives_flow) then
     allocate(aux_var_pert_dn%xmol(option%nflowspec),aux_var_pert_dn%diff(option%nflowspec))
     allocate(aux_var_pert_up%xmol(option%nflowspec),aux_var_pert_up%diff(option%nflowspec))
       ! No change in fluxe
   end select
 
-  Res(1:option%nflowspec) = fluxm(:)*option%flow_dt
-  Res(option%nflowdof) = fluxe*option%flow_dt
-
+  Res(1:option%nflowspec) = fluxm(:)
+  Res(option%nflowdof) = fluxe
+  
 end subroutine THBCFlux
 
 ! ************************************************************************** !
                           sec_dencpr, &
                           option,res_sec_heat)
 
-      r_p(iend) = r_p(iend) - res_sec_heat*option%flow_dt*volume_p(local_id)
+      r_p(iend) = r_p(iend) - res_sec_heat*volume_p(local_id)
     enddo   
     option%sec_vars_update = PETSC_FALSE
   endif
       
       if (enthalpy_flag) then
         r_p(local_id*option%nflowdof) = r_p(local_id*option%nflowdof) - hsrc1* &
-                                        option%flow_dt*volume_p(local_id)   
+                                        volume_p(local_id)   
       endif
 
       select case (source_sink%flow_condition%rate%itype)
         case(MASS_RATE_SS)
           r_p((local_id-1)*option%nflowdof + jh2o) = &
             r_p((local_id-1)*option%nflowdof + jh2o) - &
-            qsrc1*option%flow_dt*volume_p(local_id)
+            qsrc1*volume_p(local_id)
         case(HET_MASS_RATE_SS)
           qsrc1 = source_sink%flow_aux_real_var(ONE_INTEGER,iconn)/FMWH2O
           r_p((local_id-1)*option%nflowdof + jh2o) = &
             r_p((local_id-1)*option%nflowdof + jh2o) - &
-            qsrc1 *option%flow_dt*volume_p(local_id)
+            qsrc1*volume_p(local_id)
         case default
           write(string,*),source_sink%flow_condition%rate%itype
           option%io_buffer='TH mode source_sink%flow_condition%rate%itype = ' // &
       ! Update residual term associated with T
       if (qsrc1 > 0.d0) then ! injection
         r_p(local_id*option%nflowdof) = r_p(local_id*option%nflowdof) - &
-              qsrc1*aux_vars_ss(sum_connection)%h*option%flow_dt*volume_p(local_id)
+              qsrc1*aux_vars_ss(sum_connection)%h*volume_p(local_id)
       else
         ! extraction
         r_p(local_id*option%nflowdof) = r_p(local_id*option%nflowdof) - &
-              qsrc1*aux_vars(ghosted_id)%h*option%flow_dt*volume_p(local_id)
+              qsrc1*aux_vars(ghosted_id)%h*volume_p(local_id)
       endif
 
     enddo
                         option,jac_sec_heat)
                         
       Jup(option%nflowdof,2) = Jup(option%nflowdof,2) - &
-                               jac_sec_heat*option%flow_dt*volume_p(local_id)
+                               jac_sec_heat*volume_p(local_id)
     endif
                             
     ! scale by the volume of the cell
       end select
 
       if (qsrc1 > 0.d0) then ! injection
-        !dresT_dp = -qsrc1*hw_dp*option%flow_dt
-        dresT_dp = -qsrc1*aux_vars_ss(sum_connection)%dh_dp*option%flow_dt
-        ! dresT_dt = -qsrc1*hw_dt*option%flow_dt ! since tsrc1 is prescribed, there is no derivative
+        !dresT_dp = -qsrc1*hw_dp
+        dresT_dp = -qsrc1*aux_vars_ss(sum_connection)%dh_dp
+        ! dresT_dt = -qsrc1*hw_dt ! since tsrc1 is prescribed, there is no derivative
         istart = ghosted_id*option%nflowdof
         call MatSetValuesLocal(A,1,istart-1,1,istart-option%nflowdof,dresT_dp,ADD_VALUES,ierr)
       else
         ! extraction
-        dresT_dp = -qsrc1*aux_vars(ghosted_id)%dh_dp*option%flow_dt
-        dresT_dt = -qsrc1*aux_vars(ghosted_id)%dh_dt*option%flow_dt
+        dresT_dp = -qsrc1*aux_vars(ghosted_id)%dh_dp
+        dresT_dt = -qsrc1*aux_vars(ghosted_id)%dh_dt
         istart = ghosted_id*option%nflowdof
         call MatSetValuesLocal(A,1,istart-1,1,istart-option%nflowdof,dresT_dp,ADD_VALUES,ierr)
         call MatSetValuesLocal(A,1,istart-1,1,istart-1,dresT_dt,ADD_VALUES,ierr)

src/pflotran/thc.F90

 #endif                    
 #endif
 
-
+  J = J/option%flow_dt
+  
   J(option%nflowdof,:) = vol_frac_prim*J(option%nflowdof,:)
 
   if (option%numerical_derivatives_flow) then
 #endif
 #endif
 
-  Res(1:option%nflowdof-1) = mol(:)
-  Res(option%nflowdof) = vol_frac_prim*eng
+  Res(1:option%nflowdof-1) = mol(:)/option%flow_dt
+  Res(option%nflowdof) = vol_frac_prim*eng/option%flow_dt
 
 end subroutine THCAccumulation
 
   Jdn(option%nflowdof,2) = Jdn(option%nflowdof,2) + Dk*area*(-1.d0) + &
                            area*(global_aux_var_up%temp(1) - & 
                            global_aux_var_dn%temp(1))*dDk_dt_dn 
-  Jup = Jup*option%flow_dt
-  Jdn = Jdn*option%flow_dt
+
  ! note: Res is the flux contribution, for node up J = J + Jup
  !                                              dn J = J - Jdn  
 
   cond = Dk*area*(global_aux_var_up%temp(1) - global_aux_var_dn%temp(1)) 
   fluxe = fluxe + cond
 
-  Res(1:option%nflowdof-1) = fluxm(:) * option%flow_dt
-  Res(option%nflowdof) = fluxe * option%flow_dt
+  Res(1:option%nflowdof-1) = fluxm(:)
+  Res(option%nflowdof) = fluxe
  ! note: Res is the flux contribution, for node 1 R = R + Res_FL
  !                                              2 R = R - Res_FL  
  
 
   end select
 
-  Jdn = Jdn * option%flow_dt
-
   if (option%numerical_derivatives_flow) then
     allocate(aux_var_pert_dn%xmol(option%nflowspec),aux_var_pert_dn%diff(option%nflowspec))
     allocate(aux_var_pert_up%xmol(option%nflowspec),aux_var_pert_up%diff(option%nflowspec))
       ! No change in fluxe
   end select
 
-  Res(1:option%nflowspec) = fluxm(:)*option%flow_dt
-  Res(option%nflowdof) = fluxe*option%flow_dt
+  Res(1:option%nflowspec) = fluxm(:)
+  Res(option%nflowdof) = fluxe
 
 end subroutine THCBCFlux
 
                           sec_dencpr, &
                           option,res_sec_heat)
 
-      r_p(iend) = r_p(iend) - res_sec_heat*option%flow_dt*volume_p(local_id)
+      r_p(iend) = r_p(iend) - res_sec_heat*volume_p(local_id)
     enddo   
     option%sec_vars_update = PETSC_FALSE
   endif
       
       if (enthalpy_flag) then
         r_p(local_id*option%nflowdof) = r_p(local_id*option%nflowdof) - hsrc1* &
-                                          option%flow_dt*volume_p(local_id)   
+                                          volume_p(local_id)   
       endif         
 
       if (qsrc1 > 0.d0) then ! injection
-        !call wateos_noderiv(tsrc1,global_aux_vars(ghosted_id)%pres(1), &
-        !                    dw_kg,dw_mol,enth_src_h2o,option%scale,ierr)
-!           units: dw_mol [mol/dm^3]; dw_kg [kg/m^3]
-!           qqsrc = qsrc1/dw_mol ! [kmol/s (mol/dm^3 = kmol/m^3)]
+        ! units: dw_mol [mol/dm^3]; dw_kg [kg/m^3]
+        ! qqsrc = qsrc1/dw_mol ! [kmol/s (mol/dm^3 = kmol/m^3)]
         r_p((local_id-1)*option%nflowdof + jh2o) =  &
                                      r_p((local_id-1)*option%nflowdof + jh2o) &
-                                     - qsrc1*option%flow_dt*volume_p(local_id)
-        !r_p(local_id*option%nflowdof) = r_p(local_id*option%nflowdof) - qsrc1*enth_src_h2o*option%flow_dt
+                                     - qsrc1*volume_p(local_id)
         r_p(local_id*option%nflowdof) = r_p(local_id*option%nflowdof) - &
-          qsrc1*aux_vars_ss(sum_connection)%h*option%flow_dt*volume_p(local_id)
+          qsrc1*aux_vars_ss(sum_connection)%h*volume_p(local_id)
       else
         ! extraction
         r_p((local_id)*option%nflowdof+jh2o) = r_p((local_id-1)*option%nflowdof+jh2o) &
-                                               - qsrc1*option%flow_dt* &
-                                                 volume_p(local_id)
+                                               - qsrc1*volume_p(local_id)
         r_p(local_id*option%nflowdof) = r_p(local_id*option%nflowdof) - &
                                         qsrc1*aux_vars(ghosted_id)%h* &
-                                        option%flow_dt*volume_p(local_id)
+                                        volume_p(local_id)
                                         
       endif  
     
       if (csrc1 > 0.d0) then ! injection
         call printErrMsg(option,"concentration source not yet implemented in THC")
       endif
-  !  else if (qsrc1 < 0.d0) then ! withdrawal
-  !  endif
     enddo
     source_sink => source_sink%next
   enddo
                         option,jac_sec_heat)
                         
       Jup(option%nflowdof,2) = Jup(option%nflowdof,2) - &
-                               jac_sec_heat*option%flow_dt*volume_p(local_id)
+                               jac_sec_heat*volume_p(local_id)
     endif
                             
 ! scale by the volume of the cell
 
       if (associated(patch%imat)) then
         if (patch%imat(ghosted_id) <= 0) cycle
-      endif
-      
-!      if (enthalpy_flag) then
-!        r_p(local_id*option%nflowdof) = r_p(local_id*option%nflowdof) - hsrc1 * option%flow_dt   
-!      endif         
+      endif        
 
       if (qsrc1 > 0.d0) then ! injection
-        !call wateos(tsrc1,global_aux_vars(ghosted_id)%pres(1),dw_kg,dw_mol,dw_dp,dw_dt, &
-        !      enth_src_h2o,hw_dp,hw_dt,option%scale,ierr)
-!           units: dw_mol [mol/dm^3]; dw_kg [kg/m^3]
-!           qqsrc = qsrc1/dw_mol ! [kmol/s (mol/dm^3 = kmol/m^3)]
-        ! base on r_p() = r_p() - qsrc1*enth_src_h2o*option%flow_dt
-        !dresT_dp = -qsrc1*hw_dp*option%flow_dt
-        dresT_dp = -qsrc1*aux_vars_ss(sum_connection)%dh_dp*option%flow_dt
-        ! dresT_dt = -qsrc1*hw_dt*option%flow_dt ! since tsrc1 is prescribed, there is no derivative
+        ! units: dw_mol [mol/dm^3]; dw_kg [kg/m^3]
+        ! qqsrc = qsrc1/dw_mol ! [kmol/s (mol/dm^3 = kmol/m^3)]
+        ! base on r_p() = r_p() - qsrc1*enth_src_h2o
+        ! dresT_dp = -qsrc1*hw_dp
+        dresT_dp = -qsrc1*aux_vars_ss(sum_connection)%dh_dp
+        ! dresT_dt = -qsrc1*hw_dt ! since tsrc1 is prescribed, there is no derivative
         istart = ghosted_id*option%nflowdof
         call MatSetValuesLocal(A,1,istart-1,1,istart-option%nflowdof,dresT_dp,ADD_VALUES,ierr)
         ! call MatSetValuesLocal(A,1,istart-1,1,istart-1,dresT_dt,ADD_VALUES,ierr)
       else
         ! extraction
-        dresT_dp = -qsrc1*aux_vars(ghosted_id)%dh_dp*option%flow_dt
-        dresT_dt = -qsrc1*aux_vars(ghosted_id)%dh_dt*option%flow_dt
+        dresT_dp = -qsrc1*aux_vars(ghosted_id)%dh_dp
+        dresT_dt = -qsrc1*aux_vars(ghosted_id)%dh_dt
         istart = ghosted_id*option%nflowdof
         call MatSetValuesLocal(A,1,istart-1,1,istart-option%nflowdof,dresT_dp,ADD_VALUES,ierr)
         call MatSetValuesLocal(A,1,istart-1,1,istart-1,dresT_dt,ADD_VALUES,ierr)

src/pflotran/thmc.F90

                     sat_i*ddeni_dt*u_i + sat_i*den_i*dui_dt)*porXvol
 #endif
 
+  J = J/option%flow_dt 
+
   if (option%numerical_derivatives_flow) then
     allocate(thmc_aux_var_pert%xmol(option%nflowspec),thmc_aux_var_pert%diff(option%nflowspec))
     call GlobalAuxVarInit(global_aux_var_pert,option)  
   eng = eng + (sat_g*den_g*u_g + sat_i*den_i*u_i)*porXvol
 #endif
 
-  Res(1:option%nflowspec) = mol(:)
-  Res(option%nflowdof-option%nmechdof) = eng
+  Res(1:option%nflowspec) = mol(:)/option%flow_dt
+  Res(option%nflowdof-option%nmechdof) = eng/option%flow_dt
  
 
 end subroutine THMCAccumulation
                            area*(global_aux_var_up%temp(1) - & 
                            global_aux_var_dn%temp(1))*dDk_dt_dn 
                            
-  Jup = Jup*option%flow_dt
-  Jdn = Jdn*option%flow_dt
  ! note: Res is the flux contribution, for node up J = J + Jup
  !                                              dn J = J - Jdn  
 
       lambda_avg*(-dd_up*(Minv_up(1,j)*disp_vec_sum_up(1) + Minv_up(2,j)* &
       disp_vec_sum_up(2) + Minv_up(3,j)*disp_vec_sum_up(3)) + &
       dd_dn*(Minv_dn(1,j)*disp_up_minus_dn(1) + Minv_dn(2,j)* &
-      disp_up_minus_dn(2) + Minv_dn(3,j)*disp_up_minus_dn(3)))*normal(i)*area*option%flow_dt
+      disp_up_minus_dn(2) + Minv_dn(3,j)*disp_up_minus_dn(3)))*normal(i)*area
     enddo
   enddo
       
       mu_avg*(-dd_up*Minv_up(i,j)*(disp_vec_sum_up(1)*normal(1) + &
       disp_vec_sum_up(2)*normal(2) + disp_vec_sum_up(3)*normal(3)) + &
       dd_dn*Minv_dn(i,j)*(disp_up_minus_dn(1)*normal(1) + &
-      disp_up_minus_dn(2)*normal(2) + disp_up_minus_dn(3)*normal(3)))*area*option%flow_dt
+      disp_up_minus_dn(2)*normal(2) + disp_up_minus_dn(3)*normal(3)))*area
     enddo
   enddo      
       
       mu_avg*(-dd_up*disp_vec_sum_up(i)*(Minv_up(1,j)*normal(1) + &
       Minv_up(2,j)*normal(2) + Minv_up(3,j)*normal(3)) + &
       dd_dn*disp_up_minus_dn(i)*(Minv_dn(1,j)*normal(1) + &
-      Minv_dn(2,j)*normal(2) + Minv_dn(3,j)*normal(3)))*area*option%flow_dt
+      Minv_dn(2,j)*normal(2) + Minv_dn(3,j)*normal(3)))*area
     enddo
   enddo  
 
       lambda_avg*(-dd_dn*(Minv_dn(1,j)*disp_vec_sum_dn(1) + Minv_dn(2,j)* &
       disp_vec_sum_dn(2) + Minv_dn(3,j)*disp_vec_sum_dn(3)) + &
       dd_up*(Minv_up(1,j)*disp_dn_minus_up(1) + Minv_up(2,j)* &
-      disp_dn_minus_up(2) + Minv_up(3,j)*disp_dn_minus_up(3)))*normal(i)*area*option%flow_dt
+      disp_dn_minus_up(2) + Minv_up(3,j)*disp_dn_minus_up(3)))*normal(i)*area
     enddo
   enddo
       
       mu_avg*(-dd_dn*Minv_dn(i,j)*(disp_vec_sum_dn(1)*normal(1) + &
       disp_vec_sum_dn(2)*normal(2) + disp_vec_sum_dn(3)*normal(3)) + &
       dd_up*Minv_up(i,j)*(disp_dn_minus_up(1)*normal(1) + &
-      disp_dn_minus_up(2)*normal(2) + disp_dn_minus_up(3)*normal(3)))*area*option%flow_dt
+      disp_dn_minus_up(2)*normal(2) + disp_dn_minus_up(3)*normal(3)))*area
     enddo
   enddo      
       
       mu_avg*(-dd_dn*disp_vec_sum_dn(i)*(Minv_dn(1,j)*normal(1) + &
       Minv_dn(2,j)*normal(2) + Minv_dn(3,j)*normal(3)) + &
       dd_up*disp_dn_minus_up(i)*(Minv_up(1,j)*normal(1) + &
-      Minv_up(2,j)*normal(2) + Minv_up(3,j)*normal(3)))*area*option%flow_dt
+      Minv_up(2,j)*normal(2) + Minv_up(3,j)*normal(3)))*area
     enddo
   enddo  
 
   cond = Dk*area*(global_aux_var_up%temp(1) - global_aux_var_dn%temp(1)) 
   fluxe = fluxe + cond
 
-  Res(1:option%nflowspec) = fluxm(:)*option%flow_dt
-  Res(option%nflowdof-option%nmechdof) = fluxe*option%flow_dt
+  Res(1:option%nflowspec) = fluxm(:)
+  Res(option%nflowdof-option%nmechdof) = fluxe
  ! note: Res is the flux contribution, for node 1 R = R + Res_FL
  !                                              2 R = R - Res_FL  
 
                                   
   Res(option%nflowdof-option%nmechdof+1) = (stress(1,1)*unit_normal(1) + &
                                            stress(1,2)*unit_normal(2) + &
-                                           stress(1,3)*unit_normal(3))*area*option%flow_dt
+                                           stress(1,3)*unit_normal(3))*area
 
   Res(option%nflowdof-option%nmechdof+2) = (stress(2,1)*unit_normal(1) + &
                                            stress(2,2)*unit_normal(2) + &
-                                           stress(2,3)*unit_normal(3))*area*option%flow_dt
+                                           stress(2,3)*unit_normal(3))*area
 
   Res(option%nflowdof-option%nmechdof+3) = (stress(3,1)*unit_normal(1) + &
                                            stress(3,2)*unit_normal(2) + &
-                                           stress(3,3)*unit_normal(3))*area*option%flow_dt
+                                           stress(3,3)*unit_normal(3))*area
   
  
 end subroutine THMCFlux
 
   end select
 
-  Jdn = Jdn * option%flow_dt
-
   if (option%numerical_derivatives_flow) then
     allocate(aux_var_pert_dn%xmol(option%nflowspec),aux_var_pert_dn%diff(option%nflowspec))
     allocate(aux_var_pert_up%xmol(option%nflowspec),aux_var_pert_up%diff(option%nflowspec))
       ! No change in fluxe
   end select
 
-  Res(1:option%nflowspec) = fluxm(:)*option%flow_dt
-  Res(option%nflowdof-option%nmechdof) = fluxe*option%flow_dt
+  Res(1:option%nflowspec) = fluxm(:)
+  Res(option%nflowdof-option%nmechdof) = fluxe
 
 end subroutine THMCBCFlux
 
 
     r_p(iend-2) =  r_p(iend-2) + &
          thmc_parameter%rock_den(int(ithrm_loc_p(ghosted_id)))* &
-         option%gravity(1)*volume_p(local_id)*1.d-6*option%flow_dt !convert to MN
+         option%gravity(1)*volume_p(local_id)*1.d-6 !convert to MN
     r_p(iend-1) =  r_p(iend-2) + &
          thmc_parameter%rock_den(int(ithrm_loc_p(ghosted_id)))* &
-         option%gravity(2)*volume_p(local_id)*1.d-6*option%flow_dt !convert to MN
+         option%gravity(2)*volume_p(local_id)*1.d-6 !convert to MN
     r_p(iend) =  r_p(iend) + &
          thmc_parameter%rock_den(int(ithrm_loc_p(ghosted_id)))* &
-         option%gravity(3)*volume_p(local_id)*1.d-6*option%flow_dt !convert to MN
+         option%gravity(3)*volume_p(local_id)*1.d-6 !convert to MN
 
   enddo 
 
       if (enthalpy_flag) then
         r_p(local_id*(option%nflowdof-option%nmechdof)) = &
                          r_p(local_id*(option%nflowdof-option%nmechdof)) - &
-                         hsrc1*option%flow_dt*volume_p(local_id)   
+                         hsrc1*volume_p(local_id)   
       endif         
 
 !      if (qsrc1 > 0.d0) then ! injection
 !           units: dw_mol [mol/dm^3]; dw_kg [kg/m^3]
 !           qqsrc = qsrc1/dw_mol ! [kmol/s (mol/dm^3 = kmol/m^3)]
 !        r_p((local_id-1)*option%nflowdof + jh2o) = r_p((local_id-1)*option%nflowdof + jh2o) &
-!                                               - qsrc1 *option%flow_dt
-!        r_p(local_id*option%nflowdof) = r_p(local_id*option%nflowdof) - qsrc1*enth_src_h2o*option%flow_dt
+!                                               - qsrc1
+!        r_p(local_id*option%nflowdof) = r_p(local_id*option%nflowdof) - qsrc1*enth_src_h2o
 !      endif  
 !    
 !      if (csrc1 > 0.d0) then ! injection
       endif
       
 !      if (enthalpy_flag) then
-!        r_p(local_id*option%nflowdof) = r_p(local_id*option%nflowdof) - hsrc1 * option%flow_dt   
+!        r_p(local_id*option%nflowdof) = r_p(local_id*option%nflowdof) - hsrc1   
 !      endif         
 
       if (qsrc1 > 0.d0) then ! injection
               enth_src_h2o,hw_dp,hw_dt,option%scale,ierr)        
 !           units: dw_mol [mol/dm^3]; dw_kg [kg/m^3]
 !           qqsrc = qsrc1/dw_mol ! [kmol/s (mol/dm^3 = kmol/m^3)]
-        ! base on r_p() = r_p() - qsrc1*enth_src_h2o*option%flow_dt
-        dresT_dp = -qsrc1*hw_dp*option%flow_dt
-        ! dresT_dt = -qsrc1*hw_dt*option%flow_dt ! since tsrc1 is prescribed, there is no derivative
+        ! base on r_p() = r_p() - qsrc1*enth_src_h2o
+        dresT_dp = -qsrc1*hw_dp
+        ! dresT_dt = -qsrc1*hw_dt ! since tsrc1 is prescribed, there is no derivative
         istart = ghosted_id*option%nflowdof
         call MatSetValuesLocal(A,1,istart-1,1,istart-option%nflowdof,dresT_dp,ADD_VALUES,ierr)
         ! call MatSetValuesLocal(A,1,istart-1,1,istart-1,dresT_dt,ADD_VALUES,ierr)