Skip to content

Commit

Permalink
Added inviscid case, more B0 backgrounds, fixed m=2 libration forcing
Browse files Browse the repository at this point in the history
  • Loading branch information
@repepo
repepo committed Jun 24, 2022
1 parent cbcc0df commit 1316813
Show file tree
Hide file tree
Showing 27 changed files with 5,301 additions and 143 deletions.
121 changes: 80 additions & 41 deletions bin/assemble.py
100644 → 100755
View file
Original file line number Diff line number Diff line change
Expand Up @@ -158,7 +158,7 @@ def main():
row = np.arange(pos,pos+ut.N1) # start at pos+2 because of the 2 rows with bc's
col = np.zeros(ut.N1)

bdat = Ib*ut.chebco(-2, ut.N1, 3e-16, par.ricb, ut.rcmb) * ( -4*np.sqrt(5)*(par.ricb**5)/(1-par.ricb**5) ) * par.forcing_amplitude
bdat = op.Ig*ut.chebco(-2, ut.N1, 3e-16, par.ricb, ut.rcmb) * ( -4*np.sqrt(5)*(par.ricb**5)/(1-par.ricb**5) ) * par.forcing_amplitude_cmb
B = ss.csr_matrix( ( bdat, (row,col) ), shape=(ut.sizmat,1) )

np.savez('B_forced.npz', data=B.data, indices=B.indices, indptr=B.indptr, shape=B.shape)
Expand Down Expand Up @@ -416,7 +416,7 @@ def main():
print(' m=2 radial velocity forcing ')
print('--------------------------------------------')

if (par.symm == 1 and par.m==2) and (par.bci+par.bco==2):
if (par.symm == 1 and par.m==2) and (par.bci==1):

l = 2 #
L = l*(l+1)
Expand All @@ -425,9 +425,12 @@ def main():
row = np.arange(pos,pos+2)
col = np.zeros(2)

# amplitude is max libration angle
C_icb = L*par.m*par.forcing_frequency*par.forcing_amplitude_icb/par.ricb # to be checked
C_cmb = L*par.m*par.forcing_frequency*par.forcing_amplitude_cmb
#C_icb = L*par.m*par.forcing_frequency*par.forcing_amplitude_icb/par.ricb # to be checked
#C_cmb = L*par.m*par.forcing_frequency*par.forcing_amplitude_cmb

# forcing amplitude is the radial velocity amplitude
C_icb = par.forcing_amplitude_icb * par.ricb / L
C_cmb = par.forcing_amplitude_cmb * ut.rcmb / L

bdat = np.array([C_cmb, C_icb])

Expand All @@ -436,7 +439,7 @@ def main():

else:

print('This boundary flow forcing needs symm = 1 and m = 2 and bci=bco=1')
print('This boundary flow forcing needs symm = 1 and m = 2 and bci = 1')



Expand Down Expand Up @@ -894,8 +897,11 @@ def main():
# -----------------------------------------------------------
# include boundary conditions and update loc_list
bc_u_list = bc_u_spherical( l, 'section_v' )
for q in [0,1,2]:
loc_list[q]= np.concatenate( ( loc_list[q], bc_u_list[q] ) )
if (bc_u_list is None): # no bc neded in the inviscid case
pass
else:
for q in [0,1,2]:
loc_list[q]= np.concatenate( ( loc_list[q], bc_u_list[q] ) )
# -----------------------------------------------------------


Expand Down Expand Up @@ -1376,10 +1382,18 @@ def main():
def bc_u_spherical(l,loc):
'''
Spherical boundary conditions for the velocity field,
either stress-free or no-slip.
either no-penetration (for the inviscid case), stress-free or no-slip.
'''
num_rows_u = int(2 + 2*np.sign(par.ricb)) # 2 if no IC, 4 if present
num_rows_v = int(1 + 1*np.sign(par.ricb)) # 1 if no IC, 2 if present
inviscid = (par.Ek == 0) #boolean

L = l*(l+1)

if inviscid:
num_rows_u = 1
num_rows_v = 0
else:
num_rows_u = int(2 + 2*np.sign(par.ricb)) # 2 if no IC, 4 if present
num_rows_v = int(1 + 1*np.sign(par.ricb)) # 1 if no IC, 2 if present

ixu = ( par.m + 1 - ut.s )%2 #
ixv = ( par.m + ut.s )%2
Expand All @@ -1390,51 +1404,76 @@ def bc_u_spherical(l,loc):
Tbu = bv.Tb[ixu::2,:]
Tbv = bv.Tb[ixv::2,:]

if loc == 'section_u': # 2curl eqs
if loc == 'section_u': # 2curl eqs

out = ss.dok_matrix((num_rows_u, ut.N1),dtype=complex)

if par.bco == 0: # stress-free cmb
out[ 0,:] = Tbu[:,0] # P =0
out[ 1,:] = Tbu[:,2] # P''=0
if inviscid:

out[ 0,:] = Tbu[:,0] # u_r=0

else:

elif par.bco == 1: # no-slip cmb
out[ 0,:] = Tbu[:,0] # P =0
out[ 1,:] = Tbu[:,1] # P' =0
if par.bco == 0: # stress-free cmb

if par.ricb > 0 :

if par.bci == 0: # stress-free icb
out[ 2,:] = bv.Ta[:,0] # P =0
out[ 3,:] = bv.Ta[:,2] # P''=0
if par.forcing == 9: # m=2 radial forcing
out[ 0,:] = Tbu[:,0] # P = whatever we set on the B matrix
out[ 1,:] = Tbu[:,2] - (2-L)*Tbu[:,0]/ut.rcmb**2 # Nat Schaeffer's bc notes, eq. 49: P''=(2-L)*P/rcmb^2
#out[ 1,:] = Tbu[:,2] # P''=0
else:
out[ 0,:] = Tbu[:,0] # P =0
out[ 1,:] = Tbu[:,2] # P''=0

elif par.bco == 1: # no-slip cmb

if par.forcing == 9: # m=2 radial forcing
out[ 0,:] = Tbu[:,0] # P = whatever we set on the B matrix
out[ 1,:] = Tbu[:,1] + Tbu[:,0] # P' + P/r = 0
#out[ 1,:] = Tbu[:,1] # P' = 0
else:
out[ 0,:] = Tbu[:,0] # P =0
out[ 1,:] = Tbu[:,1] # P' =0

if par.ricb > 0:

elif par.bci == 1: # no-slip icb
out[ 2,:] = bv.Ta[:,0] # P =0
out[ 3,:] = bv.Ta[:,1] # P' =0
if par.bci == 0: # stress-free icb
out[ 2,:] = bv.Ta[:,0] # P =0
out[ 3,:] = bv.Ta[:,2] # P''=0

elif par.bci == 1: # no-slip icb
out[ 2,:] = bv.Ta[:,0] # P =0
out[ 3,:] = bv.Ta[:,1] # P' =0

row0 = int(ut.N1*(l-ut.m_top)/2)
col0 = int(ut.N1*(l-ut.m_top)/2)

elif loc == 'section_v': # 1curl eqs

out = ss.dok_matrix((num_rows_v, ut.N1),dtype=complex)
if inviscid:

if par.bco == 0: # stress-free cmb
out[ 0,:] = Tbv[:,1] - Tbv[:,0]/ut.rcmb # T'-(T/r)=0

elif par.bco == 1: # no-slip cmb
out[ 0,:] = Tbv[:,0] # T=0

if par.ricb > 0 :
return None

if par.bci == 0: # stress-free icb
out[ 1,:] = bv.Ta[:,1]-bv.Ta[:,0]/par.ricb # T'-(T/r)=0

elif par.bci == 1: # no-slip icb
out[ 1,:] = bv.Ta[:,0] # T=0
else:

out = ss.dok_matrix((num_rows_v, ut.N1),dtype=complex)

if par.bco == 0: # stress-free cmb
out[ 0,:] = Tbv[:,1] - Tbv[:,0]/ut.rcmb # T'-(T/r)=0

elif par.bco == 1: # no-slip cmb
out[ 0,:] = Tbv[:,0] # T=0

if par.ricb > 0 :

if par.bci == 0: # stress-free icb
out[ 1,:] = bv.Ta[:,1]-bv.Ta[:,0]/par.ricb # T'-(T/r)=0

elif par.bci == 1: # no-slip icb
out[ 1,:] = bv.Ta[:,0] # T=0

row0 = ut.n + int( ut.N1*(l-ut.m_bot)/2 )
col0 = ut.n + int( ut.N1*(l-ut.m_bot)/2 )

row0 = ut.n + int( ut.N1*(l-ut.m_bot)/2 )
col0 = ut.n + int( ut.N1*(l-ut.m_bot)/2 )

out = out.tocoo()
out2 = [out.data, out.row + row0, out.col + col0]
Expand Down
Loading

0 comments on commit 1316813

Please sign in to comment.