From cc1135469950260a467f6dd1bd5dcf3edfb1a793 Mon Sep 17 00:00:00 2001 From: Lionel Siess Date: Tue, 20 Sep 2022 23:06:20 +0200 Subject: [PATCH 1/2] (cooling) fix implicit cooling when cooling rate goes to zero --- src/main/cooling_solver.f90 | 25 +++++++++++++++---------- src/setup/setup_shock.F90 | 8 +++++--- 2 files changed, 20 insertions(+), 13 deletions(-) diff --git a/src/main/cooling_solver.f90 b/src/main/cooling_solver.f90 index 08e46e4ea..d34311c0f 100644 --- a/src/main/cooling_solver.f90 +++ b/src/main/cooling_solver.f90 @@ -142,20 +142,25 @@ 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) @@ -163,11 +168,11 @@ subroutine implicit_cooling (ui, dudt, rho, dt, mu, gamma, Tdust, K2, kappa) 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 @@ -175,10 +180,9 @@ subroutine implicit_cooling (ui, dudt, rho, dt, mu, gamma, Tdust, K2, kappa) 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) @@ -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 @@ -470,7 +475,7 @@ subroutine write_options_cooling_solver(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(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) diff --git a/src/setup/setup_shock.F90 b/src/setup/setup_shock.F90 index 8cf56b33e..4b73f568d 100644 --- a/src/setup/setup_shock.F90 +++ b/src/setup/setup_shock.F90 @@ -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 @@ -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) @@ -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) From 8cd493368561e3a84f48d0528cfde55ee759b866 Mon Sep 17 00:00:00 2001 From: Lionel Siess Date: Tue, 20 Sep 2022 23:21:22 +0200 Subject: [PATCH 2/2] (cooling_solver) fix typo --- src/main/cooling_solver.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/main/cooling_solver.f90 b/src/main/cooling_solver.f90 index d34311c0f..2cbf417a9 100644 --- a/src/main/cooling_solver.f90 +++ b/src/main/cooling_solver.f90 @@ -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)* @@ -473,7 +473,7 @@ 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 (do not modify! set by setup)',iunit) endif