Instabilities in Temperature solution-Annulus

[Prasanna2989]

Hi All,

I was running a mantle convection model in an annulus with UW2 regional Mesh and encountered with some instabilities in the Temperature solution in the Eulerian domain (Swarm is not involved). For first 2000 steps model works fine and after that it creates the instabilities as per attached images. I have tested this model with 2 different resolutions of (64x128) and (128x256). Issue exists for both cases. I herewith include the code for your reference. This code was build based on the Van Keken 2001 example in publications in UW2 regionalMesh. The model is iso-viscous. What could be the reason for this behaviour and How should I fix this?

cheers Prasanna

import underworld as uw
import glucifer
import numpy as np
from underworld import function as fn
import math
import os
import numpy as np
import h5py
import time
import mpi4py

import matplotlib.pyplot as plt
import matplotlib.pylab as pylab

comm = mpi4py.MPI.COMM_WORLD
rank = comm.Get_rank()
size = comm.Get_size()

uw.matplotlib_inline()
plt.ion()

path =os.path.join(os.path.abspath("."),"InternalHeatBenchmark/AR1/SaveData/AnnulusTestCluster2/")#On Laptop

inputPath = path
outputPath = path

if uw.mpi.rank==0:
    if not os.path.exists(outputPath):
        os.makedirs(outputPath)
uw.mpi.barrier()

# In[ ]:

# Whether or not to load the existing temperature field
loaddata = False
restart = False

# In[ ]:

#Test 32*96
annulus = uw.mesh.FeMesh_Annulus(elementRes=(128,256), 
                                  radialLengths=(0.4292,1.4292), angularExtent=(0.,360.),
                                  periodic = [False, True])

tField = uw.mesh.MeshVariable(annulus, nodeDofCount=1)
vField = uw.mesh.MeshVariable(annulus, nodeDofCount=2)
pField = uw.mesh.MeshVariable(annulus.subMesh, nodeDofCount=1)
tDotField = uw.mesh.MeshVariable(annulus, nodeDofCount=1)

mH = annulus.save(os.path.join(outputPath, "annulus.h5") )

Ra = 1e7

tempMax = 1.0
tempMin = 0.0

t_outer = tempMin
t_inner = tempMax

outer = annulus.specialSets["MaxI_VertexSet"]
inner = annulus.specialSets["MinI_VertexSet"]

boundaryNodes = outer+inner#####

# setup parameters for temperature distribution
dr = annulus.radialLengths[1] - annulus.radialLengths[0]
dT_dr = (t_outer-t_inner)/(dr)
c0 = t_inner - dT_dr*annulus.radialLengths[0]

# wavenumber for perturbation
k = 4.0#####

for ind,coord in enumerate(annulus.data):
    r = np.sqrt(coord[0]**2 + coord[1]**2)
    theta = np.arctan2(coord[1], coord[0])

    pert = 0.2 *np.sin(k*theta)

    t = r*dT_dr + c0
    tField.data[ind] = min([max([0.,t + 1.*pert]),1])

tField.data[inner.data] = t_inner
tField.data[outer.data] = t_outer

if loaddata:
    tField.load('vanKeken2001data/temp.h5',interpolate=False)#########

vBC = uw.conditions.DirichletCondition( variable=vField, indexSetsPerDof=(inner+outer, inner+outer))
tBC = uw.conditions.DirichletCondition( variable=tField, indexSetsPerDof=(inner+outer))

advDiffSLE = uw.systems.AdvectionDiffusion(tField, tDotField, vField, fn_diffusivity=1.0, conditions=tBC)

g  = 1.0* annulus._fn_unitvec_radial()
bodyForceFn = g * tField * Ra / (annulus.radialLengths[1]-annulus.radialLengths[0])

isoviscous = True
viscosity_M = 1.
delta_T_M = tempMax-tempMin

if isoviscous:
    mname = "isovisc"
else:
    mname = "stag"

if isoviscous:
    #ArrFunction = 1.
    ArrFunction = uw.function.misc.constant(1.)
else:
    E = 23.03
    ArrFunction = viscosity_M * uw.function.math.exp((E/(tField +1.))-(E/(delta_T_M +1.)))

yielding = False

strainRate_2ndInvariant = fn.tensor.second_invariant( 
                            fn.tensor.symmetric( 
                            vField.fn_gradient ))

tau0 = 1.0e5#3e5#1e5
tau1 = 1.0e7#e7#fc#7e5#1e7
tauY = tau0 + tau1*(1.-coord[1])
#tauY = tau0#Constant yielding
viscosity_Y_M = tauY/(2.*(strainRate_2ndInvariant+1.0e-18))

if yielding:
    ArrFunction = fn.exception.SafeMaths( fn.misc.min(viscosity_Y_M, ArrFunction) )

stokesSLE = uw.systems.Stokes(vField, 
                              pField, 
                              fn_viscosity=ArrFunction, 
                              fn_bodyforce=bodyForceFn,
                              conditions=vBC, 
                              _removeBCs=False)

# In[ ]:

stokesSolver = uw.systems.Solver(stokesSLE)

def update():

    cFactor = 0.5
    dt = cFactor *advDiffSLE.get_max_dt()

    advDiffSLE.integrate(dt)

    return simtime+dt, step+1

# In[ ]:

def checkpoint_function(checkpoint_number, time_Myears):

    #sH = swarm.save(os.path.join(outputPath, 'swarm-%s.h5' % checkpoint_number))

    file_prefix = os.path.join(outputPath, 'velocity-%s' % checkpoint_number)
    handle = vField.save('%s.h5' % file_prefix)
    vField.xdmf('%s.xdmf' % file_prefix, handle, 'velocity', mH, 'annulus', modeltime=time_Myears)#mH

    file_prefix = os.path.join(outputPath, 'pressure-%s' % checkpoint_number)
    handle = pField.save('%s.h5' % file_prefix)
    pField.xdmf('%s.xdmf' % file_prefix, handle, 'pressure', mH, 'annulus', modeltime=time_Myears)

    file_prefix = os.path.join(outputPath, 'temperature-%s' % checkpoint_number)
    handle = tField.save('%s.h5' % file_prefix)
    tField.xdmf('%s.xdmf' % file_prefix, handle, 'temperature', mH, 'annulus', modeltime=time_Myears)

    file_prefix = os.path.join(outputPath, 'temperatureDot-%s' % checkpoint_number)
    handle = tDotField.save('%s.h5' % file_prefix)
    tDotField.xdmf('%s.xdmf' % file_prefix, handle, 'temperatureDot', mH, 'annulus', modeltime=time_Myears)

    #file_prefix = os.path.join(outputPath, 'material-%s' % checkpoint_number)
    #handle = materialVariable.save('%s.h5' % file_prefix)
    #materialVariable.xdmf('%s.xdmf' % file_prefix, handle, 'material', sH, 'swarm', modeltime=time_Myears)

# In[ ]:

step = 0
steps_end = 3000
steps_output = 50
simtime = 0.0

# In[ ]:

time_scaling_Myr = 1.

while step < steps_end:
    stokesSolver.solve(nonLinearIterate = yielding, nonLinearMaxIterations = 15)
    simtime, step = update()
    #population_control.repopulate()
    if step % steps_output == 0:      
        #population_control.repopulate()
        checkpoint_function(step, (simtime * time_scaling_Myr))

# In[ ]:

[julesghub]

With a 1:2 ratio for the element resolution the elements geometry will be quite stretched. This can effect numerical robustness, especially running 1k+ timesteps. Consider increasing the number of "angular" elements to make the element geometry closer to a 1:1 ratio.

[Prasanna2989]

Thanks Jule. r:2 pi r (1:5) ratio models are woking fine.

Cheers