From e3247a609beac3860d9cac78d2ec5ae8161e5ded Mon Sep 17 00:00:00 2001 From: root Date: Tue, 30 Jul 2024 15:15:36 +0000 Subject: [PATCH 01/13] heatin=2 case succesfuly implemented --- .../boussinesq/shell/rtc/trace_marginal.py | 2 +- .../shell/rtc/implicit/physical_model.py | 62 ++++++++++++++++--- 2 files changed, 55 insertions(+), 9 deletions(-) diff --git a/Python/linear_stability/boussinesq/shell/rtc/trace_marginal.py b/Python/linear_stability/boussinesq/shell/rtc/trace_marginal.py index 5d53c01..387c768 100644 --- a/Python/linear_stability/boussinesq/shell/rtc/trace_marginal.py +++ b/Python/linear_stability/boussinesq/shell/rtc/trace_marginal.py @@ -149,7 +149,7 @@ Ek = (Ta*(1.0-rratio)**4)**-0.5 else: Ek = Ta**-0.5 -eq_params = {'ekman':Ek, 'prandtl':Pr, 'rayleigh':Ra, 'r_ratio':rratio, 'heating':heating} +eq_params = {'ekman':Ek, 'prandtl':Pr, 'rayleigh':Ra, 'r_ratio':rratio, 'heating':heating, 'beta':beta, 'alpha':alpha } eq_params.update(model.automatic_parameters(eq_params)) bcs = {'bcType':model.SOLVER_HAS_BC, 'velocity':bc_vel, 'temperature':bc_temp} diff --git a/Python/quicc_solver/model/boussinesq/shell/rtc/implicit/physical_model.py b/Python/quicc_solver/model/boussinesq/shell/rtc/implicit/physical_model.py index 9ff48d6..542a90a 100644 --- a/Python/quicc_solver/model/boussinesq/shell/rtc/implicit/physical_model.py +++ b/Python/quicc_solver/model/boussinesq/shell/rtc/implicit/physical_model.py @@ -22,7 +22,14 @@ def periodicity(self): def nondimensional_parameters(self): """Get the list of nondimensional parameters""" - return ["ekman", "prandtl", "rayleigh", "r_ratio", "heating"] + return ["ekman", "prandtl", "rayleigh", "r_ratio", "heating","alpha","beta"] + + # alpha = 1 yields a linear gravity in r, whereas alpha !=1 leads to a alpha*r+(1-alpha)/r^2 type gravity profile + # beta = 1 yields a linear dT/dr in r, whereas beta !=1 leads to a beta*r+(1-beta)/r^2 type dT/dr profile + # heating=0: internal heating; + # heating=1: differential heating; + # heating=2: differential heating + internal heating (proportional to r) + # heating=3: differential heating + internal heating (general form) def automatic_parameters(self, eq_params): """Extend parameters with automatically computable values""" @@ -258,8 +265,17 @@ def explicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri if field_row == ("temperature","") and field_col == ("velocity","pol"): if eq_params["heating"] == 0: mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, -bg_eff, with_sh_coeff = 'laplh', restriction = restriction) - else: + elif eq_params["heating"] == 1: mat = geo.i2(res[0], ri, ro, res[1], m, bc, -bg_eff, with_sh_coeff = 'laplh', restriction = restriction) + elif eq_params["heating"] == 2: + beta = eq_params['beta'] + # assumes length-scale is the depth of the shell + c1 = beta*bg_eff/ro # internal heating contribution + c2 = ro**2*(1-beta*bg_eff) # differential heating contribution + if beta==1: # same as for heating=0 + mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, c1, with_sh_coeff = 'laplh', restriction = restriction) + else: + mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, c1, with_sh_coeff = 'laplh', restriction = restriction) + geo.i2(res[0], ri, ro, res[1], m, bc, c2, with_sh_coeff = 'laplh', restriction = restriction) if mat is None: raise RuntimeError("Equations are not setup properly!") @@ -280,8 +296,14 @@ def nonlinear_block(self, res, eq_params, eigs, bcs, field_row, field_col, restr if field_row == ("temperature","") and field_col == field_row: if eq_params["heating"] == 0: mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, restriction = restriction) - else: + elif eq_params["heating"] == 1: mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, restriction = restriction) + elif eq_params["heating"] == 2: + beta = eq_params['beta'] + if beta==1: # same as for heating=0 + mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, restriction = restriction) + else: + mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, restriction = restriction) if mat is None: raise RuntimeError("Equations are not setup properly!") @@ -341,8 +363,17 @@ def implicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri if self.linearize: if eq_params["heating"] == 0: mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, bg_eff/Pr, with_sh_coeff = 'laplh', restriction = restriction) - else: + elif eq_params["heating"] == 1: mat = geo.i2(res[0], ri, ro, res[1], m, bc, bg_eff/Pr, with_sh_coeff = 'laplh', restriction = restriction) + elif eq_params["heating"] == 2: + beta = eq_params['beta'] + # assumes length-scale is the depth of the shell + c1 = beta*bg_eff/ro # internal heating contribution + c2 = ro**2*(1-beta*bg_eff) # differential heating contribution + if beta==1: # same as for heating=0 + mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, c1/Pr, with_sh_coeff = 'laplh', restriction = restriction) + else: + mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, c1/Pr, with_sh_coeff = 'laplh', restriction = restriction) + geo.i2(res[0], ri, ro, res[1], m, bc, c2/Pr, with_sh_coeff = 'laplh', restriction = restriction) else: mat = geo.zblk(res[0], ri, ro, res[1], m, bc) @@ -350,8 +381,14 @@ def implicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri elif field_col == ("temperature",""): if eq_params["heating"] == 0: mat = geo.i2r2lapl(res[0], ri, ro, res[1], m, bc, 1.0/(Pr*c_dt), restriction = restriction) - else: + elif eq_params["heating"] == 1: mat = geo.i2r3lapl(res[0], ri, ro, res[1], m, bc, 1.0/(Pr*c_dt), restriction = restriction) + elif eq_params["heating"] == 2: + beta = eq_params['beta'] + if beta==1: # same as for heating=0 + mat = geo.i2r2lapl(res[0], ri, ro, res[1], m, bc, 1.0/(Pr*c_dt), restriction = restriction) + else: + mat = geo.i2r3lapl(res[0], ri, ro, res[1], m, bc, 1.0/(Pr*c_dt), restriction = restriction) if mat is None: raise RuntimeError("Equations are not setup properly!") @@ -378,8 +415,14 @@ def time_block(self, res, eq_params, eigs, bcs, field_row, restriction = None): elif field_row == ("temperature",""): if eq_params["heating"] == 0: mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, restriction = restriction) - else: + elif eq_params["heating"] == 1: mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, restriction = restriction) + elif eq_params["heating"] == 2: + beta = eq_params['beta'] + if beta == 1: # same as for heating=0 + mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, restriction = restriction) + else: + mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, restriction = restriction) if mat is None: raise RuntimeError("Equations are not setup properly!") @@ -440,9 +483,12 @@ def nondimensional_factors(self, eq_params): elif eq_params['heating'] == 1: # (R_o - R_i) rescaling Ra_eff = (Ra*T/ro) - bg_eff = ro**2*rratio - + bg_eff = ro**2*rratio # Rescale Rayleigh by E^{-4/3} # Ra_eff = Ra_eff*T**(1./3.) + elif eq_params['heating'] == 2: + Ra_eff = (Ra*T/ro) + bg_eff = 1 + return (Ra_eff, bg_eff) From f4e4bafdba9606958b86527c8de681be09570c43 Mon Sep 17 00:00:00 2001 From: Stefano Maffei Date: Wed, 31 Jul 2024 19:46:06 +0200 Subject: [PATCH 02/13] fixed the rest of the python scripts --- .../boussinesq/shell/rtc/trace_marginal.py | 9 ++- .../shell/rtc/explicit/physical_model.py | 56 +++++++++++++++++-- .../shell/rtc/implicit/physical_model.py | 20 ++++--- 3 files changed, 71 insertions(+), 14 deletions(-) diff --git a/Python/linear_stability/boussinesq/shell/rtc/trace_marginal.py b/Python/linear_stability/boussinesq/shell/rtc/trace_marginal.py index 387c768..dc21194 100644 --- a/Python/linear_stability/boussinesq/shell/rtc/trace_marginal.py +++ b/Python/linear_stability/boussinesq/shell/rtc/trace_marginal.py @@ -130,6 +130,13 @@ #bc_vel = 0; bc_temp = 0; heating = 0; rratio = 0.6; Pr = 1; rescale = False #Ek = 1e-4; Rac = 27.0031; mc = 18 +# NS/NS, FT/FT, heating=2, ri/ro=0.35 +#bc_vel = 0; bc_temp = 0 +#heating = 2; alpha=0; beta=1 +#rratio = 0.35; Pr = 0.1; rescale = False +#Ek = 0.5e-3; Rac = 18.747; mc = 3 +#res = [32, 32, 0] + # Convert Ekman to Taylor number Ta = Ek**-2 @@ -179,7 +186,7 @@ def wave(m): marginal_options['save_physical'] = True marginal_options['impose_symmetry'] = False marginal_options['use_spherical_evp'] = False -marginal_options['curve_points'] = np.arange(max(1,m-3), m+10, 1) +marginal_options['curve_points'] = np.arange(max(1,m-3), m+4, 1) # Compute MarginalCurve.compute(gevp_opts, marginal_options) diff --git a/Python/quicc_solver/model/boussinesq/shell/rtc/explicit/physical_model.py b/Python/quicc_solver/model/boussinesq/shell/rtc/explicit/physical_model.py index 5d37455..6d30577 100644 --- a/Python/quicc_solver/model/boussinesq/shell/rtc/explicit/physical_model.py +++ b/Python/quicc_solver/model/boussinesq/shell/rtc/explicit/physical_model.py @@ -23,7 +23,14 @@ def periodicity(self): def nondimensional_parameters(self): """Get the list of nondimensional parameters""" - return ["ekman", "prandtl", "rayleigh", "r_ratio", "heating"] + return ["ekman", "prandtl", "rayleigh", "r_ratio", "heating","alpha","beta"] + + # alpha = 1 yields a linear gravity in r, whereas alpha !=1 leads to a alpha*r+(1-alpha)/r^2 type gravity profile + # beta = 1 yields a linear dT/dr in r, whereas beta !=1 leads to a beta*r+(1-beta)/r^2 type dT/dr profile + # heating=0: internal heating; + # heating=1: differential heating; + # heating=2: differential heating + internal heating (proportional to r) + # heating=3: differential heating + internal heating (general form) def automatic_parameters(self, eq_params): """Extend parameters with automatically computable values""" @@ -236,13 +243,26 @@ def explicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri mat = None bc = self.convert_bc(eq_params,eigs,bcs,field_row,field_col) if field_row == ("velocity","pol") and field_col == ("temperature",""): - mat = geo.i4r4(res[0], ri, ro, bc, Ra_eff) + alpha = eq_params['alpha'] + c1 = Ra_eff * alpha + c2 = Ra_eff *(1-alpha) * ro**3 + #mat = geo.i4r4(res[0], ri, ro, bc, Ra_eff) + mat = geo.i4r4(res[0], ri, ro, bc, c1) + geo.i4r1(res[0], ri, ro, bc, c2) elif field_row == ("temperature","") and field_col == ("velocity","pol"): if eq_params["heating"] == 0: mat = geo.i2r2(res[0], ri, ro, bc, -bg_eff*l*(l+1.0)) - else: + elif eq_params["heating"] == 1: mat = geo.i2(res[0], ri, ro, bc, -bg_eff*l*(l+1.0)) + elif eq_params["heating"] == 2 or eq_params["heating"] == 3: + beta = eq_params['beta'] + # assumes length-scale is the depth of the shell + c1 = beta*bg_eff/ro # internal heating contribution + c2 = ro**2*(1-beta*bg_eff) # differential heating contribution + if beta==1: # same as for heating=0 + mat = geo.i2r2(res[0], ri, ro, bc, -c1*l*(l+1.0)) + else: + mat = geo.i2r3(res[0], ri, ro, bc, -c1*l*(l+1.0)) + geo.i2(res[0], ri, ro, bc, -c2*l*(l+1.0)) if mat is None: raise RuntimeError("Equations are not setup properly!") @@ -259,8 +279,14 @@ def nonlinear_block(self, res, eq_params, eigs, bcs, field_row, field_col, restr if field_row == ("temperature","") and field_col == field_row: if eq_params["heating"] == 0: mat = geo.i2r2(res[0], ri, ro, bc) - else: + elif eq_params["heating"] == 1: mat = geo.i2r3(res[0], ri, ro, bc) + elif eq_params["heating"] == 2 or eq_params["heating"] == 3: + beta = eq_params['beta'] + if beta==1: # same as for heating=0 + mat = geo.i2r2(res[0], ri, ro, bc) + else: + mat = geo.i2r3(res[0], ri, ro, bc) if mat is None: raise RuntimeError("Equations are not setup properly!") @@ -288,8 +314,14 @@ def implicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri elif field_row == ("temperature","") and field_col == field_row: if eq_params["heating"] == 0: mat = geo.i2r2lapl(res[0], ri, ro, l, bc, 1.0/Pr) - else: + elif eq_params["heating"] == 1: mat = geo.i2r3lapl(res[0], ri, ro, l, bc, 1.0/Pr) + elif eq_params["heating"] == 2 or eq_params["heating"] == 3: + beta = eq_params['beta'] + if beta==1: # same as for heating=0 + mat = geo.i2r2lapl(res[0], ri, ro, l, bc, 1.0/Pr) + else: + mat = geo.i2r3lapl(res[0], ri, ro, l, bc, 1.0/Pr) if mat is None: raise RuntimeError("Equations are not setup properly!") @@ -315,8 +347,14 @@ def time_block(self, res, eq_params, eigs, bcs, field_row, restriction = None): elif field_row == ("temperature",""): if eq_params["heating"] == 0: mat = geo.i2r2(res[0], ri, ro, bc) - else: + elif eq_params["heating"] == 1: mat = geo.i2r3(res[0], ri, ro, bc) + elif eq_params["heating"] == 2 or eq_params["heating"] == 3: + beta = eq_params['beta'] + if beta == 1: # same as for heating=0 + mat = geo.i2r2(res[0], ri, ro, bc) + else: + mat = geo.i2r3(res[0], ri, ro, bc) if mat is None: raise RuntimeError("Equations are not setup properly!") @@ -361,5 +399,11 @@ def nondimensional_factors(self, eq_params): # (R_o - R_i) rescaling Ra_eff = (Ra*T/ro) bg_eff = ro**2*rratio + elif eq_params['heating'] == 2: + Ra_eff = (Ra*T/ro) + bg_eff = 1 + elif eq_params['heating'] == 3: + Ra_eff = (Ra*T/ro) + bg_eff = 1/(1-rratio**3) return (Ra_eff, bg_eff) diff --git a/Python/quicc_solver/model/boussinesq/shell/rtc/implicit/physical_model.py b/Python/quicc_solver/model/boussinesq/shell/rtc/implicit/physical_model.py index 542a90a..9f79230 100644 --- a/Python/quicc_solver/model/boussinesq/shell/rtc/implicit/physical_model.py +++ b/Python/quicc_solver/model/boussinesq/shell/rtc/implicit/physical_model.py @@ -267,7 +267,7 @@ def explicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, -bg_eff, with_sh_coeff = 'laplh', restriction = restriction) elif eq_params["heating"] == 1: mat = geo.i2(res[0], ri, ro, res[1], m, bc, -bg_eff, with_sh_coeff = 'laplh', restriction = restriction) - elif eq_params["heating"] == 2: + elif eq_params["heating"] == 2 or eq_params["heating"] == 3: beta = eq_params['beta'] # assumes length-scale is the depth of the shell c1 = beta*bg_eff/ro # internal heating contribution @@ -298,7 +298,7 @@ def nonlinear_block(self, res, eq_params, eigs, bcs, field_row, field_col, restr mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, restriction = restriction) elif eq_params["heating"] == 1: mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, restriction = restriction) - elif eq_params["heating"] == 2: + elif eq_params["heating"] == 2 or eq_params["heating"] == 3: beta = eq_params['beta'] if beta==1: # same as for heating=0 mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, restriction = restriction) @@ -353,7 +353,11 @@ def implicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri mat = mat + geo.i4r4lapl(res[0], ri, ro, res[1], m, bc, 1j*m*T/c_dt, with_sh_coeff = 'invlaplh', l_zero_fix = 'zero', restriction = restriction) elif field_col == ("temperature",""): - mat = geo.i4r4(res[0], ri, ro, res[1], m, bc, -Ra_eff/c_dt**2, l_zero_fix = 'zero', restriction = restriction) + alpha = eq_params['alpha'] + c1 = -Ra_eff * alpha /c_dt**2 + c2 = -Ra_eff *(1-alpha) * ro**3 /c_dt**2 + #mat = geo.i4r4(res[0], ri, ro, res[1], m, bc, -Ra_eff/c_dt**2, l_zero_fix = 'zero', restriction = restriction) + mat = geo.i4r4(res[0], ri, ro, res[1], m, bc, c1, l_zero_fix = 'zero', restriction = restriction) + geo.i4r1(res[0], ri, ro, res[1], m, bc, c2, l_zero_fix = 'zero', restriction = restriction) elif field_row == ("temperature",""): if field_col == ("velocity","tor"): @@ -365,7 +369,7 @@ def implicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, bg_eff/Pr, with_sh_coeff = 'laplh', restriction = restriction) elif eq_params["heating"] == 1: mat = geo.i2(res[0], ri, ro, res[1], m, bc, bg_eff/Pr, with_sh_coeff = 'laplh', restriction = restriction) - elif eq_params["heating"] == 2: + elif eq_params["heating"] == 2 or eq_params["heating"] == 3: beta = eq_params['beta'] # assumes length-scale is the depth of the shell c1 = beta*bg_eff/ro # internal heating contribution @@ -383,7 +387,7 @@ def implicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri mat = geo.i2r2lapl(res[0], ri, ro, res[1], m, bc, 1.0/(Pr*c_dt), restriction = restriction) elif eq_params["heating"] == 1: mat = geo.i2r3lapl(res[0], ri, ro, res[1], m, bc, 1.0/(Pr*c_dt), restriction = restriction) - elif eq_params["heating"] == 2: + elif eq_params["heating"] == 2 or eq_params["heating"] == 3: beta = eq_params['beta'] if beta==1: # same as for heating=0 mat = geo.i2r2lapl(res[0], ri, ro, res[1], m, bc, 1.0/(Pr*c_dt), restriction = restriction) @@ -417,7 +421,7 @@ def time_block(self, res, eq_params, eigs, bcs, field_row, restriction = None): mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, restriction = restriction) elif eq_params["heating"] == 1: mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, restriction = restriction) - elif eq_params["heating"] == 2: + elif eq_params["heating"] == 2 or eq_params["heating"] == 3: beta = eq_params['beta'] if beta == 1: # same as for heating=0 mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, restriction = restriction) @@ -489,6 +493,8 @@ def nondimensional_factors(self, eq_params): elif eq_params['heating'] == 2: Ra_eff = (Ra*T/ro) bg_eff = 1 - + elif eq_params['heating'] == 3: + Ra_eff = (Ra*T/ro) + bg_eff = 1/(1-rratio**3) return (Ra_eff, bg_eff) From 53f4619ee8a0de49038df433480e278ac7ddd89c Mon Sep 17 00:00:00 2001 From: Stefano Maffei Date: Wed, 31 Jul 2024 19:46:47 +0200 Subject: [PATCH 03/13] included alpha and beta in the list of parameters --- src/Model/Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp | 4 ++++ src/Model/Boussinesq/Shell/RTC/IRTCBackend.cpp | 4 ++++ src/Model/Boussinesq/Shell/RTC/IRTCModel.cpp | 2 ++ src/Model/Boussinesq/Shell/RTC/Implicit/ModelBackend.cpp | 2 ++ 4 files changed, 12 insertions(+) diff --git a/src/Model/Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp b/src/Model/Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp index 4541c18..9b017cb 100644 --- a/src/Model/Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp +++ b/src/Model/Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp @@ -33,8 +33,12 @@ #include "QuICC/PhysicalNames/Temperature.hpp" #include "QuICC/NonDimensional/Prandtl.hpp" #include "QuICC/NonDimensional/Rayleigh.hpp" +#include "QuICC/NonDimensional/Alpha.hpp" +#include "QuICC/NonDimensional/Beta.hpp" #include "QuICC/NonDimensional/Ekman.hpp" #include "QuICC/NonDimensional/Heating.hpp" +#include "QuICC/NonDimensional/Alpha.hpp" +#include "QuICC/NonDimensional/Beta.hpp" #include "QuICC/NonDimensional/RRatio.hpp" #include "QuICC/NonDimensional/Lower1d.hpp" #include "QuICC/NonDimensional/Upper1d.hpp" diff --git a/src/Model/Boussinesq/Shell/RTC/IRTCBackend.cpp b/src/Model/Boussinesq/Shell/RTC/IRTCBackend.cpp index f8645b5..62531aa 100644 --- a/src/Model/Boussinesq/Shell/RTC/IRTCBackend.cpp +++ b/src/Model/Boussinesq/Shell/RTC/IRTCBackend.cpp @@ -30,6 +30,8 @@ #include "QuICC/PhysicalNames/Temperature.hpp" #include "QuICC/NonDimensional/Prandtl.hpp" #include "QuICC/NonDimensional/Rayleigh.hpp" +#include "QuICC/NonDimensional/Alpha.hpp" +#include "QuICC/NonDimensional/Beta.hpp" #include "QuICC/NonDimensional/Ekman.hpp" #include "QuICC/NonDimensional/Heating.hpp" #include "QuICC/NonDimensional/RRatio.hpp" @@ -91,6 +93,8 @@ namespace RTC { std::vector names = { NonDimensional::Prandtl().tag(), NonDimensional::Rayleigh().tag(), + NonDimensional::Alpha().tag(), + NonDimensional::Beta().tag(), NonDimensional::Ekman().tag(), NonDimensional::Heating().tag(), NonDimensional::RRatio().tag() diff --git a/src/Model/Boussinesq/Shell/RTC/IRTCModel.cpp b/src/Model/Boussinesq/Shell/RTC/IRTCModel.cpp index 3385315..fb449da 100644 --- a/src/Model/Boussinesq/Shell/RTC/IRTCModel.cpp +++ b/src/Model/Boussinesq/Shell/RTC/IRTCModel.cpp @@ -28,6 +28,8 @@ #include "QuICC/PhysicalNames/VorticityZ.hpp" #include "QuICC/NonDimensional/Ekman.hpp" #include "QuICC/NonDimensional/Rayleigh.hpp" +#include "QuICC/NonDimensional/Alpha.hpp" +#include "QuICC/NonDimensional/Beta.hpp" #include "QuICC/NonDimensional/Prandtl.hpp" #include "QuICC/NonDimensional/RRatio.hpp" #include "QuICC/NonDimensional/Heating.hpp" diff --git a/src/Model/Boussinesq/Shell/RTC/Implicit/ModelBackend.cpp b/src/Model/Boussinesq/Shell/RTC/Implicit/ModelBackend.cpp index 3cecee3..54a0fb0 100644 --- a/src/Model/Boussinesq/Shell/RTC/Implicit/ModelBackend.cpp +++ b/src/Model/Boussinesq/Shell/RTC/Implicit/ModelBackend.cpp @@ -32,6 +32,8 @@ #include "QuICC/PhysicalNames/Temperature.hpp" #include "QuICC/NonDimensional/Prandtl.hpp" #include "QuICC/NonDimensional/Rayleigh.hpp" +#include "QuICC/NonDimensional/Alpha.hpp" +#include "QuICC/NonDimensional/Beta.hpp" #include "QuICC/NonDimensional/Ekman.hpp" #include "QuICC/NonDimensional/Heating.hpp" #include "QuICC/NonDimensional/RRatio.hpp" From f45d838b011f427f99eff8a8826eae6ecb0706d9 Mon Sep 17 00:00:00 2001 From: Stefano Maffei Date: Thu, 1 Aug 2024 20:58:22 +0200 Subject: [PATCH 04/13] general profiles implemented in the Explicit NL model --- .../shell/rtc/explicit/physical_model.py | 8 +- .../Shell/RTC/Explicit/ModelBackend.cpp | 101 ++++++++++++++++-- .../Boussinesq/Shell/RTC/IRTCBackend.cpp | 11 ++ .../Shell/RTC/Implicit/ModelBackend.cpp | 9 +- 4 files changed, 113 insertions(+), 16 deletions(-) diff --git a/Python/quicc_solver/model/boussinesq/shell/rtc/explicit/physical_model.py b/Python/quicc_solver/model/boussinesq/shell/rtc/explicit/physical_model.py index 6d30577..adb460a 100644 --- a/Python/quicc_solver/model/boussinesq/shell/rtc/explicit/physical_model.py +++ b/Python/quicc_solver/model/boussinesq/shell/rtc/explicit/physical_model.py @@ -259,10 +259,10 @@ def explicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri # assumes length-scale is the depth of the shell c1 = beta*bg_eff/ro # internal heating contribution c2 = ro**2*(1-beta*bg_eff) # differential heating contribution - if beta==1: # same as for heating=0 - mat = geo.i2r2(res[0], ri, ro, bc, -c1*l*(l+1.0)) - else: - mat = geo.i2r3(res[0], ri, ro, bc, -c1*l*(l+1.0)) + geo.i2(res[0], ri, ro, bc, -c2*l*(l+1.0)) + if beta==1: # same as for heating=0 + mat = geo.i2r2(res[0], ri, ro, bc, -c1*l*(l+1.0)) + else: + mat = geo.i2r3(res[0], ri, ro, bc, -c1*l*(l+1.0)) + geo.i2(res[0], ri, ro, bc, -c2*l*(l+1.0)) if mat is None: raise RuntimeError("Equations are not setup properly!") diff --git a/src/Model/Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp b/src/Model/Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp index 9b017cb..99eb43d 100644 --- a/src/Model/Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp +++ b/src/Model/Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp @@ -53,6 +53,7 @@ #include "QuICC/SparseSM/Chebyshev/LinearMap/I2Y3.hpp" #include "QuICC/SparseSM/Chebyshev/LinearMap/I2Y2SphLapl.hpp" #include "QuICC/SparseSM/Chebyshev/LinearMap/I2Y3SphLapl.hpp" +#include "QuICC/SparseSM/Chebyshev/LinearMap/I4Y1.hpp" #include "QuICC/SparseSM/Chebyshev/LinearMap/I4Y4.hpp" #include "QuICC/SparseSM/Chebyshev/LinearMap/I4Y4SphLapl.hpp" #include "QuICC/SparseSM/Chebyshev/LinearMap/I4Y4SphLapl2.hpp" @@ -235,17 +236,31 @@ namespace Explicit { { auto Pr = nds.find(NonDimensional::Prandtl::id())->second->value(); auto heatingMode = nds.find(NonDimensional::Heating::id())->second->value(); + auto beta = nds.find(NonDimensional::Beta::id())->second->value(); if(heatingMode == 0) { SparseSM::Chebyshev::LinearMap::I2Y2SphLapl spasm(nN, nN, ri, ro, l); decMat.real() = (1.0/Pr)*spasm.mat(); } - else + else if(heatingMode == 1) { SparseSM::Chebyshev::LinearMap::I2Y3SphLapl spasm(nN, nN, ri, ro, l); decMat.real() = (1.0/Pr)*spasm.mat(); } + else if(heatingMode == 2 || heatingMode == 3) + { + if(beta == 1) + { + SparseSM::Chebyshev::LinearMap::I2Y2SphLapl spasm(nN, nN, ri, ro, l); + decMat.real() = (1.0/Pr)*spasm.mat(); + } + else + { + SparseSM::Chebyshev::LinearMap::I2Y3SphLapl spasm(nN, nN, ri, ro, l); + decMat.real() = (1.0/Pr)*spasm.mat(); + } + } } else { @@ -284,17 +299,31 @@ namespace Explicit { else if(fieldId == std::make_pair(PhysicalNames::Temperature::id(), FieldComponents::Spectral::SCALAR)) { auto heatingMode = nds.find(NonDimensional::Heating::id())->second->value(); + auto beta = nds.find(NonDimensional::Beta::id())->second->value(); if(heatingMode == 0) { SparseSM::Chebyshev::LinearMap::I2Y2 spasm(nN, nN, ri, ro); decMat.real() = spasm.mat(); } - else + else if(heatingMode == 1) { SparseSM::Chebyshev::LinearMap::I2Y3 spasm(nN, nN, ri, ro); decMat.real() = spasm.mat(); } + else if(heatingMode == 2 || heatingMode == 3) + { + if(beta == 1) + { + SparseSM::Chebyshev::LinearMap::I2Y2 spasm(nN, nN, ri, ro); + decMat.real() = spasm.mat(); + } + else + { + SparseSM::Chebyshev::LinearMap::I2Y3 spasm(nN, nN, ri, ro); + decMat.real() = spasm.mat(); + } + } } } @@ -427,6 +456,8 @@ namespace Explicit { auto ri = nds.find(NonDimensional::Lower1d::id())->second->value(); auto ro = nds.find(NonDimensional::Upper1d::id())->second->value(); + auto alpha = nds.find(NonDimensional::Alpha::id())->second->value(); + int s = this->nBc(rowId); // Explicit linear operator if(opId == ModelOperator::ExplicitLinear::id()) @@ -435,22 +466,45 @@ namespace Explicit { colId == std::make_pair(PhysicalNames::Temperature::id(),FieldComponents::Spectral::SCALAR)) { auto Ra = effectiveRa(nds); + MHDFloat c1, c2; + c1 = Ra*alpha; + c2 = Ra*(1.-alpha)*(ro*ro*ro); if(this->useSplitEquation()) { - SparseSM::Chebyshev::LinearMap::I2Y2 spasm(nN, nN, ri, ro); - decMat.real() = Ra*spasm.mat(); + if(alpha == 1) // linear gravity + { + SparseSM::Chebyshev::LinearMap::I2Y2 spasm(nN, nN, ri, ro); + decMat.real() = Ra*spasm.mat(); + } + else //general gravity + { + // need to look better into the split equation part + //SparseSM::Chebyshev::LinearMap::I2Y2 spasm(nN, nN, ri, ro); + SparseSM::Chebyshev::LinearMap::I2Y2 spasm(nN, nN, ri, ro); + decMat.real() = Ra*spasm.mat(); + } } else { - SparseSM::Chebyshev::LinearMap::I4Y4 spasm(nN, nN, ri, ro); - decMat.real() = Ra*spasm.mat(); + if(alpha == 1) //linear gravity + { + SparseSM::Chebyshev::LinearMap::I4Y4 spasm(nN, nN, ri, ro); + decMat.real() = Ra*spasm.mat(); + } + else // general gravity + { + SparseSM::Chebyshev::LinearMap::I4Y4 i4r4(nN, nN, ri, ro); + SparseSM::Chebyshev::LinearMap::I4Y1 i4r1(nN, nN, ri, ro); + decMat.real() = c1*i4r4.mat() + c2*i4r1.mat(); + } } } else if(rowId == std::make_pair(PhysicalNames::Temperature::id(), FieldComponents::Spectral::SCALAR) && colId == std::make_pair(PhysicalNames::Velocity::id(),FieldComponents::Spectral::POL)) { auto heatingMode = nds.find(NonDimensional::Heating::id())->second->value(); + auto beta = nds.find(NonDimensional::Beta::id())->second->value(); auto bg = effectiveBg(nds); auto dl = static_cast(l); auto ll1 = dl*(dl + 1.0); @@ -460,11 +514,28 @@ namespace Explicit { SparseSM::Chebyshev::LinearMap::I2Y2 spasm(nN, nN, ri, ro); decMat.real() = -bg*ll1*spasm.mat(); } - else + else if(heatingMode == 1) { SparseSM::Chebyshev::LinearMap::I2 spasm(nN, nN, ri, ro); decMat.real() = -bg*ll1*spasm.mat(); } + else if(heatingMode == 2 || heatingMode == 3) + { + MHDFloat c1, c2; + c1 = beta*bg/ro; + c2 = ro*ro*(1-beta*bg); + if(beta == 1) + { + SparseSM::Chebyshev::LinearMap::I2Y2 spasm(nN, nN, ri, ro); + decMat.real() = -c1*ll1*spasm.mat(); + } + else + { + SparseSM::Chebyshev::LinearMap::I2Y3 i2r3(nN, nN, ri, ro); + SparseSM::Chebyshev::LinearMap::I2 i2(nN, nN, ri, ro); + decMat.real() = -c1*ll1*i2r3.mat() -c2*ll1*i2.mat(); + } + } } else { @@ -478,17 +549,31 @@ namespace Explicit { if(rowId == std::make_pair(PhysicalNames::Temperature::id(), FieldComponents::Spectral::SCALAR) && rowId == colId) { auto heatingMode = nds.find(NonDimensional::Heating::id())->second->value(); + auto beta = nds.find(NonDimensional::Beta::id())->second->value(); if(heatingMode == 0) { SparseSM::Chebyshev::LinearMap::I2Y2 spasm(nN, nN, ri, ro); decMat.real() = spasm.mat(); } - else + else if(heatingMode == 1) { SparseSM::Chebyshev::LinearMap::I2Y3 spasm(nN, nN, ri, ro); decMat.real() = spasm.mat(); } + else if(heatingMode == 2 || heatingMode == 3) + { + if(beta==1) + { + SparseSM::Chebyshev::LinearMap::I2Y2 spasm(nN, nN, ri, ro); + decMat.real() = spasm.mat(); + } + else + { + SparseSM::Chebyshev::LinearMap::I2Y3 spasm(nN, nN, ri, ro); + decMat.real() = spasm.mat(); + } + } } else { diff --git a/src/Model/Boussinesq/Shell/RTC/IRTCBackend.cpp b/src/Model/Boussinesq/Shell/RTC/IRTCBackend.cpp index 62531aa..fcbdede 100644 --- a/src/Model/Boussinesq/Shell/RTC/IRTCBackend.cpp +++ b/src/Model/Boussinesq/Shell/RTC/IRTCBackend.cpp @@ -48,6 +48,7 @@ #include "QuICC/SparseSM/Chebyshev/LinearMap/I2Y3.hpp" #include "QuICC/SparseSM/Chebyshev/LinearMap/I2Y2SphLapl.hpp" #include "QuICC/SparseSM/Chebyshev/LinearMap/I2Y3SphLapl.hpp" +#include "QuICC/SparseSM/Chebyshev/LinearMap/I4Y1.hpp" #include "QuICC/SparseSM/Chebyshev/LinearMap/I4Y4.hpp" #include "QuICC/SparseSM/Chebyshev/LinearMap/I4Y4SphLapl.hpp" #include "QuICC/SparseSM/Chebyshev/LinearMap/I4Y4SphLapl2.hpp" @@ -155,6 +156,7 @@ namespace RTC { auto ro = nds.find(NonDimensional::Upper1d::id())->second->value(); auto rratio = nds.find(NonDimensional::RRatio::id())->second->value(); auto heatingMode = nds.find(NonDimensional::Heating::id())->second->value(); + auto beta = nds.find(NonDimensional::Beta::id())->second->value(); if(ro == 1.0) { @@ -171,6 +173,15 @@ namespace RTC { { effBg = ro*ro*rratio; } + // gap with mixed heating + else if(heatingMode == 2) + { + effBg = 1.0; + } + else if(heatingMode == 2) + { + effBg = 1.0 / (1-rratio*rratio*rratio); + } return effBg; } diff --git a/src/Model/Boussinesq/Shell/RTC/Implicit/ModelBackend.cpp b/src/Model/Boussinesq/Shell/RTC/Implicit/ModelBackend.cpp index 54a0fb0..bce8b18 100644 --- a/src/Model/Boussinesq/Shell/RTC/Implicit/ModelBackend.cpp +++ b/src/Model/Boussinesq/Shell/RTC/Implicit/ModelBackend.cpp @@ -54,6 +54,7 @@ #include "QuICC/SparseSM/Chebyshev/LinearMap/I2Y3SphLapl.hpp" #include "QuICC/SparseSM/Chebyshev/LinearMap/I4Y3.hpp" #include "QuICC/SparseSM/Chebyshev/LinearMap/I4Y4D1.hpp" +#include "QuICC/SparseSM/Chebyshev/LinearMap/I4Y1.hpp" #include "QuICC/SparseSM/Chebyshev/LinearMap/I4Y4.hpp" #include "QuICC/SparseSM/Chebyshev/LinearMap/I4Y4SphLapl.hpp" #include "QuICC/SparseSM/Chebyshev/LinearMap/I4Y4SphLapl2.hpp" @@ -490,7 +491,7 @@ namespace Implicit { SparseSM::Chebyshev::LinearMap::I2Y2SphLapl spasm(nN, nN, ri, ro, l); bMat = (1.0/Pr)*spasm.mat(); } - else + else if(heatingMode == 1) { SparseSM::Chebyshev::LinearMap::I2Y3SphLapl spasm(nN, nN, ri, ro, l); bMat = (1.0/Pr)*spasm.mat(); @@ -607,7 +608,7 @@ namespace Implicit { SparseSM::Chebyshev::LinearMap::I2Y2 spasm(nN, nN, ri, ro); bMat = spasm.mat(); } - else + else if(heatingMode == 1) { SparseSM::Chebyshev::LinearMap::I2Y3 spasm(nN, nN, ri, ro); bMat = spasm.mat(); @@ -854,7 +855,7 @@ namespace Implicit { } //this->addBlock(decMat.real(), spasm.mat(), rowShift, colShift, -bg*ll1); } - else + else if(heatingMode == 1) { SparseSM::Chebyshev::LinearMap::I2 spasm(nN, nN, ri, ro); bMat = -bg*ll1*spasm.mat(); @@ -906,7 +907,7 @@ namespace Implicit { } //this->addBlock(decMat.real(), spasm.mat(), rowShift, colShift); } - else + else if(heatingMode == 1) { SparseSM::Chebyshev::LinearMap::I2Y3 spasm(nN, nN, ri, ro); bMat = spasm.mat(); From 102e5d4346e55d8b5380ea8df21b81db61fdff03 Mon Sep 17 00:00:00 2001 From: Stefano Maffei Date: Fri, 2 Aug 2024 09:53:20 +0200 Subject: [PATCH 05/13] SplitEquation versio NOT implemented --- src/Model/Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Model/Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp b/src/Model/Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp index 99eb43d..492f950 100644 --- a/src/Model/Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp +++ b/src/Model/Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp @@ -479,8 +479,8 @@ namespace Explicit { } else //general gravity { - // need to look better into the split equation part - //SparseSM::Chebyshev::LinearMap::I2Y2 spasm(nN, nN, ri, ro); + // To fully implement the splitEquation version , we need to update the nonlinear operator + // work in progress SparseSM::Chebyshev::LinearMap::I2Y2 spasm(nN, nN, ri, ro); decMat.real() = Ra*spasm.mat(); } From ed8004209673a6007d1a07ff59a939d7bc3af135 Mon Sep 17 00:00:00 2001 From: Stefano Maffei Date: Fri, 2 Aug 2024 11:16:14 +0200 Subject: [PATCH 06/13] updated the Readme --- Readme.md | 72 +++++++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 67 insertions(+), 5 deletions(-) diff --git a/Readme.md b/Readme.md index f8144bf..7dc37ec 100644 --- a/Readme.md +++ b/Readme.md @@ -1,15 +1,77 @@ # Boussinesq rotating thermal convection(RTC) in a spherical shell + The model equations are -![Model equations](https://latex.codecogs.com/svg.latex?%5Cinline%20%5Cbegin%7Balign*%7D%20%5Cleft%28%5Cpartial_t%20-%20%5CDelta%5Cright%29%5Cmathbf%7Bu%7D%20%26%20%3D%20%5Cmathbf%7Bu%7D%20%5Ctimes%20%5Cleft%28%20%5Cnabla%20%5Ctimes%20%5Cmathbf%7Bu%7D%5Cright%29%20+%20%5Cfrac%7BRa%7D%7BPr%7D%20%5CTheta%20%5Cmathbf%7Br%7D%20-%20%5Cnabla%5CPi%20%5C%5C%5B0.6cm%5D%20%5Cleft%28%5Cpartial_t%20-%20%5Cfrac%7B1%7D%7BPr%7D%5CDelta%5Cright%29%5CTheta%20%26%20%3D%20%5Cfrac%7BS%7D%7BPr%7D%20-%20%5Cmathbf%7Bu%7D%5Ccdot%5Cnabla%5CTheta%5C%5C%5B0.6cm%5D%20%5Cnabla%5Ccdot%5Cmathbf%7Bu%7D%20%26%20%3D%200%20%5Cend%7Balign*%7D) +$$ +(\partial_t -\Delta){\bf u}= {\bf u}\times(\nabla\times{\bf u})- \frac{1}{E}\hat{\bf z}\times{\bf u} -\nabla \Pi + \frac{Ra}{E} \left(\frac{d}{r_o}\right)\Theta{\bf r}\\ +$$ + +$$ +(\partial_t -\frac{1}{Pr}\Delta)\Theta = -{\bf u}\cdot\nabla\Theta + \mathcal{H}{\bf u}\cdot{\bf r} +$$ with the parameters defined as -![Nondimensional parameters](https://latex.codecogs.com/svg.latex?%5Cinline%20%5Cbegin%7Balign*%7D%20Pr%20%26%20%3D%20%5Cfrac%7B%5Cnu%7D%7B%5Ckappa%7D%5C%5C%20Ra%20%26%20%3D%20%5Cfrac%7Bg%20%5Calpha%20%5Cbeta%20r_o%5E4%7D%7B%5Cnu%5Ckappa%7D%20%5Cend%7Balign*%7D) +$$ +Pr = \frac{\nu}{\kappa};\quad E=\frac{\nu}{2\Omega d^2};\quad Ra=\frac{ \alpha g_o \mathcal{T} d}{2\nu\Omega} +$$ + +where $\nu$, $\kappa$, $\alpha$ are the kinematic viscosity, thermal conductivity, and thermal expansion coefficents; $\Omega$ is the rotation rate; $r_o$ is the outer radius of the spherical shell and $d$ is its thickness; $g_o$ is defined via the gravitational acceleration, ${\bf g} = -g_o {\bf r}/r_o$. $\mathcal{T}$ is a characteristic temperature and $\mathcal{H}$ is the effective temperature background. The spherical shell has an inner radius of $r_i$, outeer radius $r_o$ and thickness $d=r_o-r_i$. + +The model equations are non-dimensionalised using $d$ as the characteristic length-scale and the viscous time-scale $\nu/d^2$. Note that in the model equation, the factor $d/r_o$ as written in the buoyancy term, involves dimensional quantities. If the outer radius $r_o$ is already non-dimensional, that factor is simply $1/r_o$. + +The definition of both $\mathcal{T}$ and $\mathcal{H}$ depend on the mode of thermal driving (differential or internal heating, controlled by the parameter `heating`) and on the boundary conditions (fixed-temperature or fixed-flux): + +- For fixed temperature boundary conditions, the background temperature is as given in Eq. 1.6 of Dormy et al., 2004 (DOI: 10.1017/S0022112003007316). Then $\mathcal{T} = \Delta T = T_i - T_o$ (where $T$ is the temperature of the resulting reference state) and: + + $\mathcal{H}=$`bg_eff` for uniform internal heating (`heating=0`) + + $\mathcal{H}=$ `bg_eff` $/r^3$ for differential heating (`heating=1`) +where `bg_eff` is a pre-factor, function of $r_o$ and $d$. + +- **The fixed-flux case is not currently implemented cleanly**. For fixed-flux at either one of the boundaries, the relevant thermal scale is defined by a thermal gradient $\beta$, so that the characteristic temperature is $\mathcal{T}=\beta d$. The exact definition of $\beta$ (and therefore of the Rayleigh number) depends on the mode of heating, but still: + + $\mathcal{H}=$`bg_eff` for uniform internal heating (`heating=0`) + + $\mathcal{H}=$ `bg_eff` $/r^3$ for differential heating (`heating=1`) +**Currently the correct** `bg_eff` **is not implemented for fixed-flux conditions, but the code can still be used to calculate the onset of convection. The resulting Rayleigh number needs, however, to be corrected to account for this discrepancy.** + + +**At present, only uniform heating is implemented in the internally heated case.** + +The system is solved using a Toroidal/Poloidal decomposition of the velocity + +$$ +{\bf u} = \nabla\times T {\bf r} + \nabla\times\nabla\times P {\bf r} . +$$ + + +## General background profiles + +**The following functionalities are currently implemented only for the Explicit models, and not for the SplitEqaution formalism.** + +The parameters $\alpha$ (`alpha`) and $\beta$ (`beta`) can be used to set up general gravity and background temperature profiles. The following are currently implemented: + +### Gravity: +The following general gravity profile is implemented: +$$ +{\bf g} = -g_o\left[\left(\frac{r}{r_o}\right)\alpha + \left(\frac{r_o}{r}\right)^2(1-\alpha)\right]\hat{\bf r} +$$ +where, with $\alpha=1$ a linear profile is obtained, adequate for an inner core of the same density as the outer core. + +**This is the feature not implemented in the SplitEqaution formalism.** + +### Temperature: + +General background temperature profiles, $T$, containing a mix of both internal and differential heating are implemented. + +There are two options, chosed via the parameter `heating`: -The system is solved using a Toroidal/Poloidal decomposition of the velocity ![u](https://latex.codecogs.com/svg.latex?%5Cinline%20%5Cmathbf%7Bu%7D): +- `heating`=2 implements a purely linear internal heating component, resulting in: +$$ +T = T'_o\left[\left(\frac{r}{r_o}\right)\beta + \left(\frac{r_o}{r}\right)^2(1-\beta)\right] +$$ -![Toroidal/Poloidal decomposition](https://latex.codecogs.com/svg.latex?%5Cinline%20%5Cmathbf%7Bu%7D%3D%5Cmathbf%7B%5Cnabla%7D%5Ctimes%20T%20%5Cmathbf%7Br%7D%20+%20%5Cmathbf%7B%5Cnabla%7D%5Ctimes%5Cmathbf%7B%5Cnabla%7D%5Ctimes%20P%20%5Cmathbf%7Br%7D) +- `heating`=3 implements a fully general form for the internal heating component, resulting in: +$$ +\frac{dT}{dr} = T'_o\left[\frac{1}{1-\chi^3}\left(\frac{r}{r_o}\right)\beta + \left(\frac{r_o}{r}\right)^2\left(1-\frac{\beta}{1-\chi^3}\right)\right] +$$ -Note: To access the equations in codecogs online LaTeX editor replace "https://latex.codecogs.com/svg.latex?" with "https://www.codecogs.com/eqnedit.php?latex=". +In both cases, tempreature is nondimensionalised via $T'_o$, the temperature gradient at the outer boundary, $r_o$. In the above definition of the Rayleigh number, therefore: $\mathcal{T} = D T'_o$. \ No newline at end of file From ac57446811c35200061fc292560a85c9c1f6c121 Mon Sep 17 00:00:00 2001 From: Stefano Maffei Date: Fri, 2 Aug 2024 11:18:49 +0200 Subject: [PATCH 07/13] updated the Readme --- Readme.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Readme.md b/Readme.md index 7dc37ec..c45a8bb 100644 --- a/Readme.md +++ b/Readme.md @@ -74,4 +74,6 @@ $$ \frac{dT}{dr} = T'_o\left[\frac{1}{1-\chi^3}\left(\frac{r}{r_o}\right)\beta + \left(\frac{r_o}{r}\right)^2\left(1-\frac{\beta}{1-\chi^3}\right)\right] $$ +In both formulations, `heating`=1 is the pure internal heating case, and `heating`=0 is the differential heating case. + In both cases, tempreature is nondimensionalised via $T'_o$, the temperature gradient at the outer boundary, $r_o$. In the above definition of the Rayleigh number, therefore: $\mathcal{T} = D T'_o$. \ No newline at end of file From ff88c4c2515e83685517af9c418bfc1aac5cb517 Mon Sep 17 00:00:00 2001 From: Stefano Maffei Date: Sun, 4 Aug 2024 10:41:26 +0200 Subject: [PATCH 08/13] fixed a mistake in the temperaturee gradients --- .../shell/rtc/explicit/physical_model.py | 10 +++++----- .../shell/rtc/implicit/physical_model.py | 16 ++++++++-------- .../Shell/RTC/Explicit/ModelBackend.cpp | 11 +++++++---- 3 files changed, 20 insertions(+), 17 deletions(-) diff --git a/Python/quicc_solver/model/boussinesq/shell/rtc/explicit/physical_model.py b/Python/quicc_solver/model/boussinesq/shell/rtc/explicit/physical_model.py index adb460a..9b22383 100644 --- a/Python/quicc_solver/model/boussinesq/shell/rtc/explicit/physical_model.py +++ b/Python/quicc_solver/model/boussinesq/shell/rtc/explicit/physical_model.py @@ -26,7 +26,7 @@ def nondimensional_parameters(self): return ["ekman", "prandtl", "rayleigh", "r_ratio", "heating","alpha","beta"] # alpha = 1 yields a linear gravity in r, whereas alpha !=1 leads to a alpha*r+(1-alpha)/r^2 type gravity profile - # beta = 1 yields a linear dT/dr in r, whereas beta !=1 leads to a beta*r+(1-beta)/r^2 type dT/dr profile + # beta = 1 is for pure internal heating, beta =0 is for pure differential heating # heating=0: internal heating; # heating=1: differential heating; # heating=2: differential heating + internal heating (proportional to r) @@ -259,7 +259,7 @@ def explicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri # assumes length-scale is the depth of the shell c1 = beta*bg_eff/ro # internal heating contribution c2 = ro**2*(1-beta*bg_eff) # differential heating contribution - if beta==1: # same as for heating=0 + if beta==bg_eff: # same as for heating=0 mat = geo.i2r2(res[0], ri, ro, bc, -c1*l*(l+1.0)) else: mat = geo.i2r3(res[0], ri, ro, bc, -c1*l*(l+1.0)) + geo.i2(res[0], ri, ro, bc, -c2*l*(l+1.0)) @@ -283,7 +283,7 @@ def nonlinear_block(self, res, eq_params, eigs, bcs, field_row, field_col, restr mat = geo.i2r3(res[0], ri, ro, bc) elif eq_params["heating"] == 2 or eq_params["heating"] == 3: beta = eq_params['beta'] - if beta==1: # same as for heating=0 + if beta==bg_eff: # same as for heating=0 mat = geo.i2r2(res[0], ri, ro, bc) else: mat = geo.i2r3(res[0], ri, ro, bc) @@ -318,7 +318,7 @@ def implicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri mat = geo.i2r3lapl(res[0], ri, ro, l, bc, 1.0/Pr) elif eq_params["heating"] == 2 or eq_params["heating"] == 3: beta = eq_params['beta'] - if beta==1: # same as for heating=0 + if beta==bg_eff: # same as for heating=0 mat = geo.i2r2lapl(res[0], ri, ro, l, bc, 1.0/Pr) else: mat = geo.i2r3lapl(res[0], ri, ro, l, bc, 1.0/Pr) @@ -351,7 +351,7 @@ def time_block(self, res, eq_params, eigs, bcs, field_row, restriction = None): mat = geo.i2r3(res[0], ri, ro, bc) elif eq_params["heating"] == 2 or eq_params["heating"] == 3: beta = eq_params['beta'] - if beta == 1: # same as for heating=0 + if beta == bg_eff: # same as for heating=0 mat = geo.i2r2(res[0], ri, ro, bc) else: mat = geo.i2r3(res[0], ri, ro, bc) diff --git a/Python/quicc_solver/model/boussinesq/shell/rtc/implicit/physical_model.py b/Python/quicc_solver/model/boussinesq/shell/rtc/implicit/physical_model.py index 9f79230..e138ce4 100644 --- a/Python/quicc_solver/model/boussinesq/shell/rtc/implicit/physical_model.py +++ b/Python/quicc_solver/model/boussinesq/shell/rtc/implicit/physical_model.py @@ -25,7 +25,7 @@ def nondimensional_parameters(self): return ["ekman", "prandtl", "rayleigh", "r_ratio", "heating","alpha","beta"] # alpha = 1 yields a linear gravity in r, whereas alpha !=1 leads to a alpha*r+(1-alpha)/r^2 type gravity profile - # beta = 1 yields a linear dT/dr in r, whereas beta !=1 leads to a beta*r+(1-beta)/r^2 type dT/dr profile + # beta = 1 is for pure internal heating, beta=0 is for pure differential heating # heating=0: internal heating; # heating=1: differential heating; # heating=2: differential heating + internal heating (proportional to r) @@ -270,9 +270,9 @@ def explicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri elif eq_params["heating"] == 2 or eq_params["heating"] == 3: beta = eq_params['beta'] # assumes length-scale is the depth of the shell - c1 = beta*bg_eff/ro # internal heating contribution - c2 = ro**2*(1-beta*bg_eff) # differential heating contribution - if beta==1: # same as for heating=0 + c1 = -beta*bg_eff/ro # internal heating contribution + c2 = -ro**2*(1-beta*bg_eff) # differential heating contribution + if beta==bg_eff: # same as for heating=0 mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, c1, with_sh_coeff = 'laplh', restriction = restriction) else: mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, c1, with_sh_coeff = 'laplh', restriction = restriction) + geo.i2(res[0], ri, ro, res[1], m, bc, c2, with_sh_coeff = 'laplh', restriction = restriction) @@ -300,7 +300,7 @@ def nonlinear_block(self, res, eq_params, eigs, bcs, field_row, field_col, restr mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, restriction = restriction) elif eq_params["heating"] == 2 or eq_params["heating"] == 3: beta = eq_params['beta'] - if beta==1: # same as for heating=0 + if beta==bg_eff: # same as for heating=0 mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, restriction = restriction) else: mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, restriction = restriction) @@ -374,7 +374,7 @@ def implicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri # assumes length-scale is the depth of the shell c1 = beta*bg_eff/ro # internal heating contribution c2 = ro**2*(1-beta*bg_eff) # differential heating contribution - if beta==1: # same as for heating=0 + if beta==bg_eff: # same as for heating=0 mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, c1/Pr, with_sh_coeff = 'laplh', restriction = restriction) else: mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, c1/Pr, with_sh_coeff = 'laplh', restriction = restriction) + geo.i2(res[0], ri, ro, res[1], m, bc, c2/Pr, with_sh_coeff = 'laplh', restriction = restriction) @@ -389,7 +389,7 @@ def implicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri mat = geo.i2r3lapl(res[0], ri, ro, res[1], m, bc, 1.0/(Pr*c_dt), restriction = restriction) elif eq_params["heating"] == 2 or eq_params["heating"] == 3: beta = eq_params['beta'] - if beta==1: # same as for heating=0 + if beta==bg_eff: # same as for heating=0 mat = geo.i2r2lapl(res[0], ri, ro, res[1], m, bc, 1.0/(Pr*c_dt), restriction = restriction) else: mat = geo.i2r3lapl(res[0], ri, ro, res[1], m, bc, 1.0/(Pr*c_dt), restriction = restriction) @@ -423,7 +423,7 @@ def time_block(self, res, eq_params, eigs, bcs, field_row, restriction = None): mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, restriction = restriction) elif eq_params["heating"] == 2 or eq_params["heating"] == 3: beta = eq_params['beta'] - if beta == 1: # same as for heating=0 + if beta == bg_eff: # same as for heating=0 mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, restriction = restriction) else: mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, restriction = restriction) diff --git a/src/Model/Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp b/src/Model/Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp index 492f950..f419671 100644 --- a/src/Model/Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp +++ b/src/Model/Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp @@ -237,6 +237,7 @@ namespace Explicit { auto Pr = nds.find(NonDimensional::Prandtl::id())->second->value(); auto heatingMode = nds.find(NonDimensional::Heating::id())->second->value(); auto beta = nds.find(NonDimensional::Beta::id())->second->value(); + auto bg = effectiveBg(nds); if(heatingMode == 0) { @@ -250,7 +251,7 @@ namespace Explicit { } else if(heatingMode == 2 || heatingMode == 3) { - if(beta == 1) + if(beta == bg) { SparseSM::Chebyshev::LinearMap::I2Y2SphLapl spasm(nN, nN, ri, ro, l); decMat.real() = (1.0/Pr)*spasm.mat(); @@ -300,6 +301,7 @@ namespace Explicit { { auto heatingMode = nds.find(NonDimensional::Heating::id())->second->value(); auto beta = nds.find(NonDimensional::Beta::id())->second->value(); + auto bg = effectiveBg(nds); if(heatingMode == 0) { @@ -313,7 +315,7 @@ namespace Explicit { } else if(heatingMode == 2 || heatingMode == 3) { - if(beta == 1) + if(beta == bg) { SparseSM::Chebyshev::LinearMap::I2Y2 spasm(nN, nN, ri, ro); decMat.real() = spasm.mat(); @@ -524,7 +526,7 @@ namespace Explicit { MHDFloat c1, c2; c1 = beta*bg/ro; c2 = ro*ro*(1-beta*bg); - if(beta == 1) + if(beta == bg) { SparseSM::Chebyshev::LinearMap::I2Y2 spasm(nN, nN, ri, ro); decMat.real() = -c1*ll1*spasm.mat(); @@ -550,6 +552,7 @@ namespace Explicit { { auto heatingMode = nds.find(NonDimensional::Heating::id())->second->value(); auto beta = nds.find(NonDimensional::Beta::id())->second->value(); + auto bg = effectiveBg(nds); if(heatingMode == 0) { @@ -563,7 +566,7 @@ namespace Explicit { } else if(heatingMode == 2 || heatingMode == 3) { - if(beta==1) + if(beta==bg) { SparseSM::Chebyshev::LinearMap::I2Y2 spasm(nN, nN, ri, ro); decMat.real() = spasm.mat(); From d389e1a886dc4c64f14f38dd57b3043317f44b17 Mon Sep 17 00:00:00 2001 From: Stefano Maffei Date: Sun, 4 Aug 2024 10:53:22 +0200 Subject: [PATCH 09/13] fixed a mistake in the temperaturee gradients --- .../boussinesq/shell/rtc/explicit/physical_model.py | 8 ++++---- .../boussinesq/shell/rtc/implicit/physical_model.py | 10 +++++----- .../Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp | 8 ++++---- 3 files changed, 13 insertions(+), 13 deletions(-) diff --git a/Python/quicc_solver/model/boussinesq/shell/rtc/explicit/physical_model.py b/Python/quicc_solver/model/boussinesq/shell/rtc/explicit/physical_model.py index 9b22383..82329f9 100644 --- a/Python/quicc_solver/model/boussinesq/shell/rtc/explicit/physical_model.py +++ b/Python/quicc_solver/model/boussinesq/shell/rtc/explicit/physical_model.py @@ -259,7 +259,7 @@ def explicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri # assumes length-scale is the depth of the shell c1 = beta*bg_eff/ro # internal heating contribution c2 = ro**2*(1-beta*bg_eff) # differential heating contribution - if beta==bg_eff: # same as for heating=0 + if beta==1/bg_eff: # same as for heating=0 mat = geo.i2r2(res[0], ri, ro, bc, -c1*l*(l+1.0)) else: mat = geo.i2r3(res[0], ri, ro, bc, -c1*l*(l+1.0)) + geo.i2(res[0], ri, ro, bc, -c2*l*(l+1.0)) @@ -283,7 +283,7 @@ def nonlinear_block(self, res, eq_params, eigs, bcs, field_row, field_col, restr mat = geo.i2r3(res[0], ri, ro, bc) elif eq_params["heating"] == 2 or eq_params["heating"] == 3: beta = eq_params['beta'] - if beta==bg_eff: # same as for heating=0 + if beta==1/bg_eff: # same as for heating=0 mat = geo.i2r2(res[0], ri, ro, bc) else: mat = geo.i2r3(res[0], ri, ro, bc) @@ -318,7 +318,7 @@ def implicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri mat = geo.i2r3lapl(res[0], ri, ro, l, bc, 1.0/Pr) elif eq_params["heating"] == 2 or eq_params["heating"] == 3: beta = eq_params['beta'] - if beta==bg_eff: # same as for heating=0 + if beta==1/bg_eff: # same as for heating=0 mat = geo.i2r2lapl(res[0], ri, ro, l, bc, 1.0/Pr) else: mat = geo.i2r3lapl(res[0], ri, ro, l, bc, 1.0/Pr) @@ -351,7 +351,7 @@ def time_block(self, res, eq_params, eigs, bcs, field_row, restriction = None): mat = geo.i2r3(res[0], ri, ro, bc) elif eq_params["heating"] == 2 or eq_params["heating"] == 3: beta = eq_params['beta'] - if beta == bg_eff: # same as for heating=0 + if beta == 1/bg_eff: # same as for heating=0 mat = geo.i2r2(res[0], ri, ro, bc) else: mat = geo.i2r3(res[0], ri, ro, bc) diff --git a/Python/quicc_solver/model/boussinesq/shell/rtc/implicit/physical_model.py b/Python/quicc_solver/model/boussinesq/shell/rtc/implicit/physical_model.py index e138ce4..23bcf81 100644 --- a/Python/quicc_solver/model/boussinesq/shell/rtc/implicit/physical_model.py +++ b/Python/quicc_solver/model/boussinesq/shell/rtc/implicit/physical_model.py @@ -272,7 +272,7 @@ def explicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri # assumes length-scale is the depth of the shell c1 = -beta*bg_eff/ro # internal heating contribution c2 = -ro**2*(1-beta*bg_eff) # differential heating contribution - if beta==bg_eff: # same as for heating=0 + if beta==1/bg_eff: # same as for heating=0 mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, c1, with_sh_coeff = 'laplh', restriction = restriction) else: mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, c1, with_sh_coeff = 'laplh', restriction = restriction) + geo.i2(res[0], ri, ro, res[1], m, bc, c2, with_sh_coeff = 'laplh', restriction = restriction) @@ -300,7 +300,7 @@ def nonlinear_block(self, res, eq_params, eigs, bcs, field_row, field_col, restr mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, restriction = restriction) elif eq_params["heating"] == 2 or eq_params["heating"] == 3: beta = eq_params['beta'] - if beta==bg_eff: # same as for heating=0 + if beta==1/bg_eff: # same as for heating=0 mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, restriction = restriction) else: mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, restriction = restriction) @@ -374,7 +374,7 @@ def implicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri # assumes length-scale is the depth of the shell c1 = beta*bg_eff/ro # internal heating contribution c2 = ro**2*(1-beta*bg_eff) # differential heating contribution - if beta==bg_eff: # same as for heating=0 + if beta==1/bg_eff: # same as for heating=0 mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, c1/Pr, with_sh_coeff = 'laplh', restriction = restriction) else: mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, c1/Pr, with_sh_coeff = 'laplh', restriction = restriction) + geo.i2(res[0], ri, ro, res[1], m, bc, c2/Pr, with_sh_coeff = 'laplh', restriction = restriction) @@ -389,7 +389,7 @@ def implicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri mat = geo.i2r3lapl(res[0], ri, ro, res[1], m, bc, 1.0/(Pr*c_dt), restriction = restriction) elif eq_params["heating"] == 2 or eq_params["heating"] == 3: beta = eq_params['beta'] - if beta==bg_eff: # same as for heating=0 + if beta==1/bg_eff: # same as for heating=0 mat = geo.i2r2lapl(res[0], ri, ro, res[1], m, bc, 1.0/(Pr*c_dt), restriction = restriction) else: mat = geo.i2r3lapl(res[0], ri, ro, res[1], m, bc, 1.0/(Pr*c_dt), restriction = restriction) @@ -423,7 +423,7 @@ def time_block(self, res, eq_params, eigs, bcs, field_row, restriction = None): mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, restriction = restriction) elif eq_params["heating"] == 2 or eq_params["heating"] == 3: beta = eq_params['beta'] - if beta == bg_eff: # same as for heating=0 + if beta == 1/bg_eff: # same as for heating=0 mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, restriction = restriction) else: mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, restriction = restriction) diff --git a/src/Model/Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp b/src/Model/Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp index f419671..183403c 100644 --- a/src/Model/Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp +++ b/src/Model/Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp @@ -251,7 +251,7 @@ namespace Explicit { } else if(heatingMode == 2 || heatingMode == 3) { - if(beta == bg) + if(beta == 1/bg) { SparseSM::Chebyshev::LinearMap::I2Y2SphLapl spasm(nN, nN, ri, ro, l); decMat.real() = (1.0/Pr)*spasm.mat(); @@ -315,7 +315,7 @@ namespace Explicit { } else if(heatingMode == 2 || heatingMode == 3) { - if(beta == bg) + if(beta == 1/bg) { SparseSM::Chebyshev::LinearMap::I2Y2 spasm(nN, nN, ri, ro); decMat.real() = spasm.mat(); @@ -526,7 +526,7 @@ namespace Explicit { MHDFloat c1, c2; c1 = beta*bg/ro; c2 = ro*ro*(1-beta*bg); - if(beta == bg) + if(beta == 1/bg) { SparseSM::Chebyshev::LinearMap::I2Y2 spasm(nN, nN, ri, ro); decMat.real() = -c1*ll1*spasm.mat(); @@ -566,7 +566,7 @@ namespace Explicit { } else if(heatingMode == 2 || heatingMode == 3) { - if(beta==bg) + if(beta==1/bg) { SparseSM::Chebyshev::LinearMap::I2Y2 spasm(nN, nN, ri, ro); decMat.real() = spasm.mat(); From 93f0c68151fce80a1f3198edddc47972f26fa58d Mon Sep 17 00:00:00 2001 From: Stefano Maffei Date: Sun, 4 Aug 2024 23:34:09 +0200 Subject: [PATCH 10/13] distinct treatement for purely differential and purely internal heating --- .../shell/rtc/explicit/physical_model.py | 4 +++- .../shell/rtc/implicit/physical_model.py | 14 +++++++++++--- .../Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp | 5 +++++ 3 files changed, 19 insertions(+), 4 deletions(-) diff --git a/Python/quicc_solver/model/boussinesq/shell/rtc/explicit/physical_model.py b/Python/quicc_solver/model/boussinesq/shell/rtc/explicit/physical_model.py index 82329f9..b08166a 100644 --- a/Python/quicc_solver/model/boussinesq/shell/rtc/explicit/physical_model.py +++ b/Python/quicc_solver/model/boussinesq/shell/rtc/explicit/physical_model.py @@ -259,8 +259,10 @@ def explicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri # assumes length-scale is the depth of the shell c1 = beta*bg_eff/ro # internal heating contribution c2 = ro**2*(1-beta*bg_eff) # differential heating contribution - if beta==1/bg_eff: # same as for heating=0 + if beta==1/bg_eff: # similar to heating=0 mat = geo.i2r2(res[0], ri, ro, bc, -c1*l*(l+1.0)) + elif beta==0: # similar to heating=1 + mat = geo.i2(res[0], ri, ro, bc, -c2*l*(l+1.0)) else: mat = geo.i2r3(res[0], ri, ro, bc, -c1*l*(l+1.0)) + geo.i2(res[0], ri, ro, bc, -c2*l*(l+1.0)) diff --git a/Python/quicc_solver/model/boussinesq/shell/rtc/implicit/physical_model.py b/Python/quicc_solver/model/boussinesq/shell/rtc/implicit/physical_model.py index 23bcf81..a988271 100644 --- a/Python/quicc_solver/model/boussinesq/shell/rtc/implicit/physical_model.py +++ b/Python/quicc_solver/model/boussinesq/shell/rtc/implicit/physical_model.py @@ -272,8 +272,10 @@ def explicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri # assumes length-scale is the depth of the shell c1 = -beta*bg_eff/ro # internal heating contribution c2 = -ro**2*(1-beta*bg_eff) # differential heating contribution - if beta==1/bg_eff: # same as for heating=0 + if beta==1/bg_eff: # similar to heating=0 mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, c1, with_sh_coeff = 'laplh', restriction = restriction) + elif beta==0: # similar to heating==1 + mat = geo.i2(res[0], ri, ro, res[1], m, bc, c2, with_sh_coeff = 'laplh', restriction = restriction) else: mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, c1, with_sh_coeff = 'laplh', restriction = restriction) + geo.i2(res[0], ri, ro, res[1], m, bc, c2, with_sh_coeff = 'laplh', restriction = restriction) @@ -291,6 +293,8 @@ def nonlinear_block(self, res, eq_params, eigs, bcs, field_row, field_col, restr ri, ro = (self.automatic_parameters(eq_params)['lower1d'], self.automatic_parameters(eq_params)['upper1d']) + Ra_eff, bg_eff = self.nondimensional_factors(eq_params) + mat = None bc = self.convert_bc(eq_params,eigs,bcs,field_row,field_col) if field_row == ("temperature","") and field_col == field_row: @@ -300,7 +304,7 @@ def nonlinear_block(self, res, eq_params, eigs, bcs, field_row, field_col, restr mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, restriction = restriction) elif eq_params["heating"] == 2 or eq_params["heating"] == 3: beta = eq_params['beta'] - if beta==1/bg_eff: # same as for heating=0 + if beta==1/bg_eff: # similar to heating=0 mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, restriction = restriction) else: mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, restriction = restriction) @@ -374,8 +378,10 @@ def implicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri # assumes length-scale is the depth of the shell c1 = beta*bg_eff/ro # internal heating contribution c2 = ro**2*(1-beta*bg_eff) # differential heating contribution - if beta==1/bg_eff: # same as for heating=0 + if beta==1/bg_eff: # similar to heating=0 mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, c1/Pr, with_sh_coeff = 'laplh', restriction = restriction) + elif beta==0: # similar to heating=1 + mat = geo.i2(res[0], ri, ro, res[1], m, bc, c2/Pr, with_sh_coeff = 'laplh', restriction = restriction) else: mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, c1/Pr, with_sh_coeff = 'laplh', restriction = restriction) + geo.i2(res[0], ri, ro, res[1], m, bc, c2/Pr, with_sh_coeff = 'laplh', restriction = restriction) @@ -408,6 +414,8 @@ def time_block(self, res, eq_params, eigs, bcs, field_row, restriction = None): ri, ro = (self.automatic_parameters(eq_params)['lower1d'], self.automatic_parameters(eq_params)['upper1d']) + Ra_eff, bg_eff = self.nondimensional_factors(eq_params) + mat = None bc = self.convert_bc(eq_params,eigs,bcs,field_row,field_row) if field_row == ("velocity","tor"): diff --git a/src/Model/Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp b/src/Model/Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp index 183403c..ed01894 100644 --- a/src/Model/Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp +++ b/src/Model/Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp @@ -531,6 +531,11 @@ namespace Explicit { SparseSM::Chebyshev::LinearMap::I2Y2 spasm(nN, nN, ri, ro); decMat.real() = -c1*ll1*spasm.mat(); } + else if(beta == 0) + { + SparseSM::Chebyshev::LinearMap::I2 spasm(nN, nN, ri, ro); + decMat.real() = -c2*ll1*spasm.mat(); + } else { SparseSM::Chebyshev::LinearMap::I2Y3 i2r3(nN, nN, ri, ro); From e74c0d2e91246118c05320c7ed608931ad145617 Mon Sep 17 00:00:00 2001 From: Stefano Maffei Date: Mon, 5 Aug 2024 11:25:22 +0200 Subject: [PATCH 11/13] fixed an error --- .../model/boussinesq/shell/rtc/explicit/physical_model.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/Python/quicc_solver/model/boussinesq/shell/rtc/explicit/physical_model.py b/Python/quicc_solver/model/boussinesq/shell/rtc/explicit/physical_model.py index b08166a..886ee05 100644 --- a/Python/quicc_solver/model/boussinesq/shell/rtc/explicit/physical_model.py +++ b/Python/quicc_solver/model/boussinesq/shell/rtc/explicit/physical_model.py @@ -276,6 +276,8 @@ def nonlinear_block(self, res, eq_params, eigs, bcs, field_row, field_col, restr ri, ro = (self.automatic_parameters(eq_params)['lower1d'], self.automatic_parameters(eq_params)['upper1d']) + Ra_eff, bg_eff = self.nondimensional_factors(eq_params) + mat = None bc = self.convert_bc(eq_params,eigs,bcs,field_row,field_col) if field_row == ("temperature","") and field_col == field_row: @@ -305,6 +307,8 @@ def implicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri ri, ro = (self.automatic_parameters(eq_params)['lower1d'], self.automatic_parameters(eq_params)['upper1d']) + Ra_eff, bg_eff = self.nondimensional_factors(eq_params) + mat = None bc = self.convert_bc(eq_params,eigs,bcs,field_row,field_col) if field_row == ("velocity","tor") and field_col == field_row: @@ -338,6 +342,8 @@ def time_block(self, res, eq_params, eigs, bcs, field_row, restriction = None): ri, ro = (self.automatic_parameters(eq_params)['lower1d'], self.automatic_parameters(eq_params)['upper1d']) + Ra_eff, bg_eff = self.nondimensional_factors(eq_params) + mat = None bc = self.convert_bc(eq_params,eigs,bcs,field_row,field_row) if field_row == ("velocity","tor"): From 24239861a73e63cf9ce81349c09c2cb2e8526e9d Mon Sep 17 00:00:00 2001 From: Stefano Maffei Date: Thu, 22 Aug 2024 21:35:49 +0200 Subject: [PATCH 12/13] Fixed documentation, Fixed Pr dependency and fixed-flux conditions in the python implicit model --- .../shell/rtc/implicit/physical_model.py | 14 +++++++++----- Readme.md | 4 ++++ 2 files changed, 13 insertions(+), 5 deletions(-) diff --git a/Python/quicc_solver/model/boussinesq/shell/rtc/implicit/physical_model.py b/Python/quicc_solver/model/boussinesq/shell/rtc/implicit/physical_model.py index a988271..d2654c0 100644 --- a/Python/quicc_solver/model/boussinesq/shell/rtc/implicit/physical_model.py +++ b/Python/quicc_solver/model/boussinesq/shell/rtc/implicit/physical_model.py @@ -199,12 +199,16 @@ def convert_bc(self, eq_params, eigs, bcs, field_row, field_col): bc = {0:-22, 'rt':0} elif field_col == ("velocity","pol"): bc = {0:-41, 'rt':0} + elif field_col == ("temperature",""): + bc = {0:-21, 'rt':0} else: if field_row == ("velocity","tor") and field_col == ("velocity","tor"): bc = {0:22} elif field_row == ("velocity","pol") and field_col == ("velocity","pol"): bc = {0:41} + elif field_row == ("temperature","") and field_col == ("temperature",""): + bc = {0:21} # Set LHS galerkin restriction if self.use_galerkin: @@ -370,20 +374,20 @@ def implicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri elif field_col == ("velocity","pol"): if self.linearize: if eq_params["heating"] == 0: - mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, bg_eff/Pr, with_sh_coeff = 'laplh', restriction = restriction) + mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, bg_eff, with_sh_coeff = 'laplh', restriction = restriction) elif eq_params["heating"] == 1: - mat = geo.i2(res[0], ri, ro, res[1], m, bc, bg_eff/Pr, with_sh_coeff = 'laplh', restriction = restriction) + mat = geo.i2(res[0], ri, ro, res[1], m, bc, bg_eff, with_sh_coeff = 'laplh', restriction = restriction) elif eq_params["heating"] == 2 or eq_params["heating"] == 3: beta = eq_params['beta'] # assumes length-scale is the depth of the shell c1 = beta*bg_eff/ro # internal heating contribution c2 = ro**2*(1-beta*bg_eff) # differential heating contribution if beta==1/bg_eff: # similar to heating=0 - mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, c1/Pr, with_sh_coeff = 'laplh', restriction = restriction) + mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, c1, with_sh_coeff = 'laplh', restriction = restriction) elif beta==0: # similar to heating=1 - mat = geo.i2(res[0], ri, ro, res[1], m, bc, c2/Pr, with_sh_coeff = 'laplh', restriction = restriction) + mat = geo.i2(res[0], ri, ro, res[1], m, bc, c2, with_sh_coeff = 'laplh', restriction = restriction) else: - mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, c1/Pr, with_sh_coeff = 'laplh', restriction = restriction) + geo.i2(res[0], ri, ro, res[1], m, bc, c2/Pr, with_sh_coeff = 'laplh', restriction = restriction) + mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, c1, with_sh_coeff = 'laplh', restriction = restriction) + geo.i2(res[0], ri, ro, res[1], m, bc, c2, with_sh_coeff = 'laplh', restriction = restriction) else: mat = geo.zblk(res[0], ri, ro, res[1], m, bc) diff --git a/Readme.md b/Readme.md index c45a8bb..91a942d 100644 --- a/Readme.md +++ b/Readme.md @@ -51,9 +51,11 @@ The parameters $\alpha$ (`alpha`) and $\beta$ (`beta`) can be used to set up gen ### Gravity: The following general gravity profile is implemented: + $$ {\bf g} = -g_o\left[\left(\frac{r}{r_o}\right)\alpha + \left(\frac{r_o}{r}\right)^2(1-\alpha)\right]\hat{\bf r} $$ + where, with $\alpha=1$ a linear profile is obtained, adequate for an inner core of the same density as the outer core. **This is the feature not implemented in the SplitEqaution formalism.** @@ -65,11 +67,13 @@ General background temperature profiles, $T$, containing a mix of both internal There are two options, chosed via the parameter `heating`: - `heating`=2 implements a purely linear internal heating component, resulting in: + $$ T = T'_o\left[\left(\frac{r}{r_o}\right)\beta + \left(\frac{r_o}{r}\right)^2(1-\beta)\right] $$ - `heating`=3 implements a fully general form for the internal heating component, resulting in: + $$ \frac{dT}{dr} = T'_o\left[\frac{1}{1-\chi^3}\left(\frac{r}{r_o}\right)\beta + \left(\frac{r_o}{r}\right)^2\left(1-\frac{\beta}{1-\chi^3}\right)\right] $$ From 4e887e02817c05de7a126d84391ac6b6632e7a09 Mon Sep 17 00:00:00 2001 From: Stefano Maffei Date: Wed, 28 Aug 2024 11:53:09 +0200 Subject: [PATCH 13/13] fixed a bug in IRTCBackend.cpp --- src/Model/Boussinesq/Shell/RTC/IRTCBackend.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Model/Boussinesq/Shell/RTC/IRTCBackend.cpp b/src/Model/Boussinesq/Shell/RTC/IRTCBackend.cpp index fcbdede..02fec8b 100644 --- a/src/Model/Boussinesq/Shell/RTC/IRTCBackend.cpp +++ b/src/Model/Boussinesq/Shell/RTC/IRTCBackend.cpp @@ -178,7 +178,7 @@ namespace RTC { { effBg = 1.0; } - else if(heatingMode == 2) + else if(heatingMode == 3) { effBg = 1.0 / (1-rratio*rratio*rratio); }