Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

fix bug in implicit cooling when du/dt goes to zero #328

Merged
merged 3 commits into from
Sep 21, 2022
Merged

fix bug in implicit cooling when du/dt goes to zero #328

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
cc11354
(cooling) fix implicit cooling when cooling rate goes to zero
lsiess Sep 20, 2022
55f8608
Merge remote-tracking branch 'upstream/master'
lsiess Sep 20, 2022
8cd4933
(cooling_solver) fix typo
lsiess Sep 20, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
29 changes: 17 additions & 12 deletions src/main/cooling_solver.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ module cooling_solver
! - bowen_Cprime : *radiative cooling rate (g.s/cm³)*
! - dust_collision : *dust collision (1=on/0=off)*
! - excitation_HI : *cooling via electron excitation of HI (1=on/0=off)*
! - lambda_shock : *Cooling rate parmaeter for analytic shock solution*
! - lambda_shock : *Cooling rate parameter for analytic shock solution*
! - relax_bowen : *Bowen (diffusive) relaxation (1=on/0=off)*
! - relax_stefan : *radiative relaxation (1=on/0=off)*
! - shock_problem : *piecewise formulation for analytic shock solution (1=on/0=off)*
Expand Down Expand Up @@ -142,43 +142,47 @@ subroutine implicit_cooling (ui, dudt, rho, dt, mu, gamma, Tdust, K2, kappa)
real, intent(in) :: Tdust, K2, kappa
real, intent(out) :: dudt

real, parameter :: tol = 1.d-5, Tmin = 1.
real, parameter :: tol = 1.d-6, Tmin = 1.
integer, parameter :: iter_max = 40
real :: u,Q,dlnQ_dlnT,T_on_u,Qi,f0,fi,fmid,T,Ti,T0,dx,Tmid
real :: u,Q,dlnQ_dlnT,T_on_u,Qi,f0,fi,fmid,T,T0,dx,Tmid
integer :: iter

u = ui
T_on_u = (gamma-1.)*mu*unit_ergg/Rg
T = ui*T_on_u
call calc_cooling_rate(Q,dlnQ_dlnT, rho, T, Tdust, mu, gamma, K2, kappa)
!cooling negligible, return
if (abs(Q) < tiny(0.)) then
dudt = 0.
return
endif
T0 = T
f0 = -Q*dt*T_on_u
fi = f0
iter = 0
!define bisection interval f(T) = T^(n+1)-T^n-Q*dt*T_on_u
!define bisection interval for function f(T) = T^(n+1)-T^n-Q*dt*T_on_u
do while (((f0 > 0. .and. fi > 0.) .or. (f0 < 0. .and. fi < 0.)) .and. iter < iter_max)
Tmid = max(T+Q*dt*T_on_u,Tmin)
call calc_cooling_rate(Qi,dlnQ_dlnT, rho, Tmid, Tdust, mu, gamma, K2, kappa)
fi = Tmid-T0-Qi*dt*T_on_u
T = Tmid
iter = iter+1
enddo
!Temperature is between T0 and Tmid
if (iter > iter_max) stop '[implicit_cooling] cannot bracket cooling function'
iter = 0
if (Ti > T0) then
if (Tmid > T0) then
T = T0
fmid = fi
else
if (Townsend_test) then
!special fix for Townsend benchmark
T = max(Tcap,Tmid)
else
T = Tmid
endif
fmid = f0
endif
dx = abs(Ti-T0)
do while (abs(dx/T0) > tol .and. iter < iter_max)
dx = abs(Tmid-T0)
do while (dx/T0 > tol .and. iter < iter_max)
dx = dx*.5
Tmid = T+dx
call calc_cooling_rate(Qi,dlnQ_dlnT, rho, Tmid, Tdust, mu, gamma, K2, kappa)
Expand All @@ -189,7 +193,8 @@ subroutine implicit_cooling (ui, dudt, rho, dt, mu, gamma, Tdust, K2, kappa)
else
if (fmid <= 0.) T = Tmid
endif
iter = iter + 1
iter = iter + 1
!print *,iter,fmid,T,Tmid
enddo
u = Tmid/T_on_u
dudt =(u-ui)/dt
Expand Down Expand Up @@ -468,9 +473,9 @@ subroutine write_options_cooling_solver(iunit)
call write_inopt(dust_collision,'dust_collision','dust collision (1=on/0=off)',iunit)
call write_inopt(shock_problem,'shock_problem','piecewise formulation for analytic shock solution (1=on/0=off)',iunit)
if (shock_problem == 1) then
call write_inopt(lambda_shock_cgs,'lambda_shock','Cooling rate parmaeter for analytic shock solution',iunit)
call write_inopt(lambda_shock_cgs,'lambda_shock','Cooling rate parameter for analytic shock solution',iunit)
call write_inopt(T1_factor,'T1_factor','factor by which T0 is increased (T1= T1_factor*T0)',iunit)
call write_inopt(T0_value,'T0','temperature to cool towards',iunit)
call write_inopt(T0_value,'T0','temperature to cool towards (do not modify! set by setup)',iunit)
endif
call write_inopt(bowen_Cprime,'bowen_Cprime','radiative cooling rate (g.s/cm³)',iunit)

Expand Down
8 changes: 5 additions & 3 deletions src/setup/setup_shock.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -104,7 +104,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,
use setshock, only:set_shock,adjust_shock_boundaries,fsmooth
use radiation_utils, only:radiation_and_gas_temperature_equal
use eos_idealplusrad,only:get_idealgasplusrad_tempfrompres,get_idealplusrad_enfromtemp
use eos, only:gmw,temperature_coef,init_eos
use eos, only:temperature_coef,init_eos
use cooling, only:T0_value
#ifdef NONIDEALMHD
use nicil, only:eta_constant,eta_const_type,icnstsemi
Expand Down Expand Up @@ -739,7 +739,8 @@ subroutine write_setupfile(filename,iprint,numstates,gamma,polyk,dtg)
call write_inopt(tmax,'tmax','maximum runtime',lu,ierr1)
call write_inopt(dtmax,'dtmax','time between dumps',lu,ierr1)
call write_inopt(ieos,'ieos','equation of state option',lu,ierr1)
if (do_radiation) call write_inopt(gmw,'gmw','mean molecular weight',lu,ierr1)
call write_inopt(gmw,'gmw','mean molecular weight',lu,ierr1)
!if (do_radiation) call write_inopt(gmw,'gmw','mean molecular weight',lu,ierr1)
#ifdef NONIDEALMHD
call write_inopt(use_ohm,'use_ohm','include Ohmic resistivity',lu,ierr1)
call write_inopt(use_hall,'use_hall','include the Hall effect',lu,ierr1)
Expand Down Expand Up @@ -799,7 +800,8 @@ subroutine read_setupfile(filename,iprint,numstates,gamma,polyk,dtg,ierr)
call read_inopt(tmax,'tmax',db,errcount=nerr)
call read_inopt(dtmax,'dtmax',db,errcount=nerr)
call read_inopt(ieos,'ieos',db,errcount=nerr)
if (do_radiation) call read_inopt(gmw,'gmw',db,errcount=nerr)
call read_inopt(gmw,'gmw',db,errcount=nerr)
!if (do_radiation) call read_inopt(gmw,'gmw',db,errcount=nerr)
#ifdef NONIDEALMHD
call read_inopt(use_ohm,'use_ohm',db,errcount=nerr)
call read_inopt(use_hall,'use_hall',db,errcount=nerr)
Expand Down