Skip to content

[WIP] Temperature gravity general #27

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

Open
smaffei wants to merge 13 commits into
base: main
Choose a base branch
from
Open

[WIP] Temperature gravity general #27

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
e3247a6
heatin=2 case succesfuly implemented
smaffei Jul 30, 2024
f4e4baf
fixed the rest of the python scripts
Jul 31, 2024
53f4619
included alpha and beta in the list of parameters
Jul 31, 2024
f45d838
general profiles implemented in the Explicit NL model
Aug 1, 2024
102e5d4
SplitEquation versio NOT implemented
Aug 2, 2024
ed80042
updated the Readme
Aug 2, 2024
ac57446
updated the Readme
Aug 2, 2024
ff88c4c
fixed a mistake in the temperaturee gradients
Aug 4, 2024
d389e1a
fixed a mistake in the temperaturee gradients
Aug 4, 2024
93f0c68
distinct treatement for purely differential and purely internal heating
Aug 4, 2024
e74c0d2
fixed an error
Aug 5, 2024
2423986
Fixed documentation, Fixed Pr dependency and fixed-flux conditions in…
Aug 22, 2024
4e887e0
fixed a bug in IRTCBackend.cpp
Aug 28, 2024
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
11 changes: 9 additions & 2 deletions Python/linear_stability/boussinesq/shell/rtc/trace_marginal.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -149,7 +156,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}

Expand Down Expand Up @@ -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)
64 changes: 58 additions & 6 deletions Python/quicc_solver/model/boussinesq/shell/rtc/explicit/physical_model.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 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)
# heating=3: differential heating + internal heating (general form)

def automatic_parameters(self, eq_params):
"""Extend parameters with automatically computable values"""
Expand Down Expand Up @@ -236,13 +243,28 @@ 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/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))

if mat is None:
raise RuntimeError("Equations are not setup properly!")
Expand All @@ -254,13 +276,21 @@ 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:
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/bg_eff: # 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!")
Expand All @@ -277,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:
Expand All @@ -288,8 +320,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/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)

if mat is None:
raise RuntimeError("Equations are not setup properly!")
Expand All @@ -304,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"):
Expand All @@ -315,8 +355,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/bg_eff: # 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!")
Expand Down Expand Up @@ -361,5 +407,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)
86 changes: 75 additions & 11 deletions Python/quicc_solver/model/boussinesq/shell/rtc/implicit/physical_model.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 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)
# heating=3: differential heating + internal heating (general form)

def automatic_parameters(self, eq_params):
"""Extend parameters with automatically computable values"""
Expand Down Expand Up @@ -192,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:
Expand Down Expand Up @@ -258,8 +269,19 @@ 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 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, 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)

if mat is None:
raise RuntimeError("Equations are not setup properly!")
Expand All @@ -275,13 +297,21 @@ 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:
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 or eq_params["heating"] == 3:
beta = eq_params['beta']
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)

if mat is None:
raise RuntimeError("Equations are not setup properly!")
Expand Down Expand Up @@ -331,7 +361,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"):
Expand All @@ -340,18 +374,35 @@ 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)
else:
mat = geo.i2(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, 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, 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)

else:
mat = geo.zblk(res[0], ri, ro, res[1], m, bc)

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 or eq_params["heating"] == 3:
beta = eq_params['beta']
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)

if mat is None:
raise RuntimeError("Equations are not setup properly!")
Expand All @@ -367,6 +418,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"):
Expand All @@ -378,8 +431,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 or eq_params["heating"] == 3:
beta = eq_params['beta']
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)

if mat is None:
raise RuntimeError("Equations are not setup properly!")
Expand Down Expand Up @@ -440,9 +499,14 @@ 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
elif eq_params['heating'] == 3:
Ra_eff = (Ra*T/ro)
bg_eff = 1/(1-rratio**3)

return (Ra_eff, bg_eff)
Loading