From 857bb4aa2254bc04fa1dd9505bda3ee321ddbed2 Mon Sep 17 00:00:00 2001 From: Kushal Nagaraj Date: Thu, 24 Apr 2025 10:38:19 +0200 Subject: [PATCH 1/3] Update preCalculations.f90 Commenting out line 419, to adapt to cold boundary conditions on both boundaries. --- src/preCalculations.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/preCalculations.f90 b/src/preCalculations.f90 index c7c53ac1..78ce9058 100644 --- a/src/preCalculations.f90 +++ b/src/preCalculations.f90 @@ -416,7 +416,7 @@ subroutine preCalc(tscheme) if ( ktops == 1 .and. kbots == 1 .and. ( tops(0,0) == 0.0_cp .and. bots(0,0) == 0.0_cp ) ) then ! Fixed entropy tops(0,0)=0.0_cp - bots(0,0)=sq4pi + !bots(0,0)=sq4pi else if ( ktops == 3 .and. kbots == 3 ) then ! Fixed temperature contrast From ff57ff277dc82483cae93247f56945fb8f063837 Mon Sep 17 00:00:00 2001 From: Yann Gaillard Date: Thu, 31 Jul 2025 06:38:49 +0200 Subject: [PATCH 2/3] Add the volumetric heating term The volumetric heating term is enabled using the boolean l_vol_heat. This term scales with the inverse of r^2 and is meant for the use using buoyancy with the Rayleigh number. --- src/Namelists.f90 | 8 +++++--- src/get_nl.f90 | 6 +++++- src/get_td.f90 | 10 +++++----- src/logic.f90 | 1 + src/rIter.f90 | 4 ++-- 5 files changed, 18 insertions(+), 11 deletions(-) diff --git a/src/Namelists.f90 b/src/Namelists.f90 index 84a7c3a7..b8f386ff 100644 --- a/src/Namelists.f90 +++ b/src/Namelists.f90 @@ -85,9 +85,9 @@ subroutine readNamelists(tscheme) & strat,polind,DissNb,g0,g1,g2,r_cut_model,thickStrat, & & epsS,slopeStrat,rStrat,ampStrat, & & r_LCR,nVarDiff,nVarVisc,difExp,nVarEps,interior_model, & - & nVarEntropyGrad,l_isothermal,ktopp,po,prec_angle, & - & dilution_fac,stef,tmelt,phaseDiffFac,penaltyFac, & - & epsPhase,ktopphi,kbotphi,ampForce + & nVarEntropyGrad,l_isothermal,l_vol_heat,ktopp,po, & + & prec_angle,dilution_fac,stef,tmelt,phaseDiffFac, & + & penaltyFac,epsPhase,ktopphi,kbotphi,ampForce namelist/B_external/ & & rrMP,amp_imp,expo_imp,bmax_imp,n_imp,l_imp, & @@ -990,6 +990,7 @@ subroutine writeNamelists(n_out) write(n_out,'('' radratio ='',ES14.6,'','')') radratio write(n_out,'('' l_isothermal ='',l3,'','')') l_isothermal + write(n_out,'('' l_vol_heat ='',l3,'','')') l_vol_heat write(n_out,'('' phaseDiffFac ='',ES14.6,'','')') phaseDiffFac write(n_out,'('' epsPhase ='',ES14.6,'','')') epsPhase write(n_out,'('' penaltyFac ='',ES14.6,'','')') penaltyFac @@ -1361,6 +1362,7 @@ subroutine defaultNamelists l_non_rot =.false. ! No Coriolis force ! l_anel =.false. ! Anelastic stuff ! l_isothermal =.false. ! Isothermal = 0 Grünesein ! + l_vol_heat =.false. ! No volemtric heating interior_model="None" ! Name of the interior model !---- Run time and number of threads: diff --git a/src/get_nl.f90 b/src/get_nl.f90 index ae522878..2e0d2d68 100644 --- a/src/get_nl.f90 +++ b/src/get_nl.f90 @@ -30,7 +30,7 @@ module grid_space_arrays_mod use parallel_mod, only: get_openmp_blocks use constants, only: two, third, one use logic, only: l_conv_nl, l_heat_nl, l_mag_nl, l_anel, l_mag_LF, l_adv_curl, & - & l_chemical_conv, l_precession, l_centrifuge, l_phase_field, l_ehd_dep, l_ehd_die + & l_chemical_conv, l_precession, l_centrifuge, l_phase_field, l_ehd_dep, l_ehd_die, l_vol_heat implicit none @@ -460,6 +460,10 @@ subroutine get_nl(this, time, nR, nBc, lRmsCalc) this%heatTerms(:,nPhi)= & & opr * rae/rat * radratio**2/(1.0D0-radratio)**4 * or4(nR) end if + + if ( l_vol_heat ) then + this%heatTerms(:,nPhi)= opr + end if end do !$omp end parallel diff --git a/src/get_td.f90 b/src/get_td.f90 index d4a53242..955b1638 100644 --- a/src/get_td.f90 +++ b/src/get_td.f90 @@ -11,7 +11,7 @@ module nonlinear_lm_mod use truncation, only: lm_max, lm_maxMag, m_min use logic, only : l_anel, l_conv_nl, l_corr, l_heat, l_anelastic_liquid, & & l_mag_nl, l_mag_kin, l_mag_LF, l_conv, l_mag, & - & l_chemical_conv, l_single_matrix, l_double_curl, l_ehd_die + & l_chemical_conv, l_single_matrix, l_double_curl, l_ehd_die, l_vol_heat use radial_functions, only: r, or2, or1, beta, epscProf, or4, temp0, orho1, l_R use physical_parameters, only: CorFac, epsc, n_r_LCR, epscXi use blocking, only: lm2l, lm2lmA, lm2lmS @@ -60,7 +60,7 @@ subroutine initialize(this,lmP_max) allocate( this%VxBrLM(lmP_max), this%VxBtLM(lmP_max), this%VxBpLM(lmP_max)) bytes_allocated = bytes_allocated + 6*lmP_max*SIZEOF_DEF_COMPLEX - if ( l_anel .or. l_ehd_die ) then + if ( l_anel .or. l_ehd_die .or. l_vol_heat ) then allocate( this%heatTermsLM(lmP_max) ) bytes_allocated = bytes_allocated+lmP_max*SIZEOF_DEF_COMPLEX end if @@ -120,7 +120,7 @@ subroutine set_zero(this) this%VStLM(lm)=zero this%VSpLM(lm)=zero end if - if ( l_anel .or. l_ehd_die ) this%heatTermsLM(lm)=zero + if ( l_anel .or. l_ehd_die .or. l_vol_heat ) this%heatTermsLM(lm)=zero if ( l_chemical_conv ) then this%VXitLM(lm)=zero this%VXipLM(lm)=zero @@ -478,7 +478,7 @@ subroutine get_dsdt(this, nR, nBc, dsdt, dVSrLM) if ( nBc == 0 ) then dsdt(1)=epsc*epscProf(nR)!+opr/epsS*divKtemp0(nR) - if ( l_anel .or. l_ehd_die ) then + if ( l_anel .or. l_ehd_die .or. l_vol_heat ) then if ( l_anelastic_liquid ) then dsdt(1)=dsdt(1)+temp0(nR)*this%heatTermsLM(1) else @@ -486,7 +486,7 @@ subroutine get_dsdt(this, nR, nBc, dsdt, dVSrLM) end if end if - if ( l_anel .or. l_ehd_die ) then + if ( l_anel .or. l_ehd_die .or. l_vol_heat ) then if ( l_anelastic_liquid ) then !$omp parallel do do lm=lm_min,lm_max diff --git a/src/logic.f90 b/src/logic.f90 index fc295ae7..c571c6d4 100644 --- a/src/logic.f90 +++ b/src/logic.f90 @@ -74,6 +74,7 @@ module logic logical :: lVerbose ! Switch for detailed information about run progress logical :: l_ehd_dep ! Switch for dilectrophoretic force logical :: l_ehd_die ! Switch for dielectric heating + logical :: l_vol_heat ! Switch for volumetric heating logical :: l_PressGraph ! Store pressure in graphic files diff --git a/src/rIter.f90 b/src/rIter.f90 index 92c8a49d..2e14aa71 100644 --- a/src/rIter.f90 +++ b/src/rIter.f90 @@ -21,7 +21,7 @@ module rIter_mod & l_precession, l_centrifuge, l_adv_curl, & & l_double_curl, l_parallel_solve, l_single_matrix,& & l_temperature_diff, l_RMS, l_phase_field, & - & l_onset, l_DTrMagSpec, l_ehd_dep, l_ehd_die + & l_onset, l_DTrMagSpec, l_ehd_dep, l_ehd_die, l_vol_heat use radial_data, only: n_r_cmb, n_r_icb, nRstart, nRstop, nRstartMag, & & nRstopMag use radial_functions, only: or2, orho1, l_R @@ -683,7 +683,7 @@ subroutine transform_to_lm_space(this, nR, lRmsCalc, dVSrLM, dVXirLM, dphidt) call spat_to_qst(this%gsa%VSr, this%gsa%VSt, this%gsa%VSp, & & dVSrLM, this%nl_lm%VStLM, this%nl_lm%VSpLM, l_R(nR)) - if ( l_anel .or. l_ehd_die ) call scal_to_SH(this%gsa%heatTerms, this%nl_lm%heatTermsLM, & + if ( l_anel .or. l_ehd_die .or. l_vol_heat ) call scal_to_SH(this%gsa%heatTerms, this%nl_lm%heatTermsLM, & & l_R(nR)) end if if ( l_chemical_conv ) then From bfc4990dcdf512cbab893525aff2cdfb9ee439dc Mon Sep 17 00:00:00 2001 From: Kushal Date: Wed, 8 Oct 2025 13:28:41 +0200 Subject: [PATCH 3/3] Correction in the volumetric heating get_nl.f90 chnage in equation (energy), Pr replace with 1 --- src/get_nl.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/get_nl.f90 b/src/get_nl.f90 index 2e0d2d68..d3b030c7 100644 --- a/src/get_nl.f90 +++ b/src/get_nl.f90 @@ -462,7 +462,7 @@ subroutine get_nl(this, time, nR, nBc, lRmsCalc) end if if ( l_vol_heat ) then - this%heatTerms(:,nPhi)= opr + this%heatTerms(:,nPhi)= one end if end do !$omp end parallel