Merge pull request #428 from StephaneMichoulier/implicitforce

2-fluid drag implicit scheme
danieljprice · Jun 14, 2023 · b637b43 · b637b43
2 parents f486a4b + 482e4e0
commit b637b43
Show file tree

Hide file tree

Showing 10 changed files with 277 additions and 106 deletions.
diff --git a/src/main/config.F90 b/src/main/config.F90
@@ -130,7 +130,7 @@ module dim
                                    radensumden
 
  ! fsum
- integer, parameter :: fsumvars = 20 ! Number of scalars in fsum
+ integer, parameter :: fsumvars = 23 ! Number of scalars in fsum
  integer, parameter :: fsumarrs = 5  ! Number of arrays  in fsum
  integer, parameter :: maxfsum  = fsumvars + &                  ! Total number of values
                                   fsumarrs*(maxdusttypes-1) + &
@@ -139,7 +139,7 @@ module dim
 ! xpartveci
  integer, parameter :: maxxpartvecidens = 14 + radenxpartvetden
 
- integer, parameter :: maxxpartvecvars = 57 ! Number of scalars in xpartvec
+ integer, parameter :: maxxpartvecvars = 61 ! Number of scalars in xpartvec
  integer, parameter :: maxxpartvecarrs = 2  ! Number of arrays in xpartvec
  integer, parameter :: maxxpartvecGR   = 33 ! Number of GR values in xpartvec (1 for dens, 16 for gcov, 16 for gcon)
  integer, parameter :: maxxpartveciforce = maxxpartvecvars + &              ! Total number of values

diff --git a/src/main/deriv.F90 b/src/main/deriv.F90
@@ -47,7 +47,7 @@ subroutine derivs(icall,npart,nactive,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,&
  use ptmass,         only:ipart_rhomax,ptmass_calc_enclosed_mass,ptmass_boundary_crossing
  use externalforces, only:externalforce
  use part,           only:dustgasprop,dvdx,Bxyz,set_boundaries_to_active,&
-                          nptmass,xyzmh_ptmass,sinks_have_heating,dust_temp,VrelVf
+                          nptmass,xyzmh_ptmass,sinks_have_heating,dust_temp,VrelVf,fxyz_drag
 #ifdef IND_TIMESTEPS
  use timestep_ind,   only:nbinmax
 #else
@@ -184,7 +184,7 @@ subroutine derivs(icall,npart,nactive,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,&
  stressmax = 0.
  if (sinks_have_heating(nptmass,xyzmh_ptmass)) call ptmass_calc_enclosed_mass(nptmass,npart,xyzh)
  call force(icall,npart,xyzh,vxyzu,fxyzu,divcurlv,divcurlB,Bevol,dBevol,&
-            rad,drad,radprop,dustprop,dustgasprop,dustfrac,ddustevol,&
+            rad,drad,radprop,dustprop,dustgasprop,dustfrac,ddustevol,fext,fxyz_drag,&
             ipart_rhomax,dt,stressmax,eos_vars,dens,metrics)
  call do_timing('force',tlast,tcpulast)
 

diff --git a/src/main/dust.F90 b/src/main/dust.F90
@@ -40,6 +40,7 @@ module dust
  integer, public  :: icut_backreaction    = 0
  integer, public  :: irecon               = 1
  logical, public  :: ilimitdustflux       = .false. ! to limit spurious dust generation in outer disc
+ logical, public  :: drag_implicit        = .false.  ! use implicit scheme for 2-fluids drag forces
 
  public :: get_ts
  public :: init_drag
@@ -321,6 +322,7 @@ subroutine write_options_dust(iunit)
     call write_inopt(ilimitdustflux,'ilimitdustflux','limit the dust flux using Ballabio et al. (2018)',iunit)
  else
     call write_inopt(irecon,'irecon','use reconstruction in gas/dust drag (-1=off,0=no slope limiter,1=van leer MC)',iunit)
+    !call write_inopt(drag_implicit,'drag_implicit','gas/dust drag implicit scheme (!!! Works only with IND_TIMESTEPS=no !!!)',iunit)
  endif
 
  call write_inopt(icut_backreaction,'icut_backreaction','cut the drag on the gas phase (0=no, 1=yes)',iunit)
@@ -375,6 +377,8 @@ subroutine read_options_dust(name,valstring,imatch,igotall,ierr)
     !--no longer a compulsory parameter
  case('irecon')
     read(valstring,*,iostat=ierr) irecon
+ !case('drag_implicit')
+ !   read(valstring,*,iostat=ierr) drag_implicit
  case default
     imatch = .false.
  end select
@@ -399,7 +403,7 @@ subroutine read_options_dust(name,valstring,imatch,igotall,ierr)
  case(0,1)
     ineed(iKcode) = 0
  case(2,3)
-    ineed(iKcode) = 1
+    ineed(iKcode) = 0 !1
  case default
     call fatal('read_dust_infile_options','Invalid option',var='idrag',ival=idrag)
  end select