<jupyter_start><jupyter_text>Benchmark: Moresi and Solomatov 1995=======Stagnant lid convection model------**This example covers the concepts of:**1. Changing rheologies, specifically temperature dependent viscosity function.2. Stokes solver options.3. Saving and loading FE variables.**Keywords:** particle swarms, Stokes system, advective diffusive systems**References**<NAME>., and <NAME> (1995), Numerical investigation of 2d convection with extremely large viscosity variations, Phys. Fluids, 7, 2154–2162.<!--1. <NAME> and <NAME>, Mantle convection with a brittle lithosphere: thoughts on the global tectonic styles of the Earth and Venus. Geophys. J. Int. (1998) 133 (3): 669-682.  http://gji.oxfordjournals.org/content/133/3/669.short--><jupyter_code>%matplotlib inline
import matplotlib.pyplot as pyplot
import matplotlib.pylab as pylab
import numpy as np
import underworld as uw
import math
from underworld import function as fn
import glucifer
import time<jupyter_output><empty_output><jupyter_text>Setup parameters-----<jupyter_code>dim = 2
boxLength = 1.0
boxHeight = 1.0
tempMin = 0.0
tempMax = 1.0
# Set the resolution.
res = 32
# Set the Rayleigh number.
Ra=1.e6<jupyter_output><empty_output><jupyter_text>**Input/output paths**Set input and output file directory paths. For this example the input directory contains near steady state snapshots of the velocity, pressure and temperature fields. It also constains a summary of the Nusselt and $v_{RMS}$ values against time from the simulation used to make these snapshots.<jupyter_code>inputPath  = 'MandS_Input/'
outputPath = 'MandS_Output/'
# Make output directory if necessary
import os
if not os.path.exists(outputPath):
    os.makedirs(outputPath)<jupyter_output><empty_output><jupyter_text>Create mesh and finite element variables------In this case the mesh type used is different. For more information on the different mesh types see the user guide.<jupyter_code>mesh = uw.mesh.FeMesh_Cartesian( elementType = ("Q2/dPc1"), 
                                 elementRes  = (res, res), 
                                 minCoord    = (0., 0.), 
                                 maxCoord    = (boxLength, boxHeight) )
 
velocityField    = uw.mesh.MeshVariable( mesh=mesh,    nodeDofCount=dim )
pressureField    = uw.mesh.MeshVariable( mesh=mesh.subMesh,    nodeDofCount=1 )
temperatureField = uw.mesh.MeshVariable( mesh=mesh, nodeDofCount=1 )
tempDotField     = uw.mesh.MeshVariable( mesh=mesh, nodeDofCount=1 )<jupyter_output><empty_output><jupyter_text>Set initial conditions and boundary conditions----------**Initial and boundary conditions**Either set by perturbation function or load data from file.<jupyter_code>LoadFromFile = False<jupyter_output><empty_output><jupyter_text>**If loading from file**Read 32\*32 resolution data for $P$, $v$ and $T$ fields as well as existing summary statistics data. These are converted into lists so that the main time loop below will append with new values.<jupyter_code>if(LoadFromFile == True):
# set up mesh for 32*32 data file
    mesh32 = uw.mesh.FeMesh_Cartesian( elementType = ("Q2/dPc1"), 
                                               elementRes = (32, 32), 
                                                 minCoord = (0., 0.), 
                                                 maxCoord = (boxLength, boxHeight)  )
    temperatureField32 = uw.mesh.MeshVariable( mesh=mesh32,         nodeDofCount=1 )
    pressureField32    = uw.mesh.MeshVariable( mesh=mesh32.subMesh, nodeDofCount=1 )
    velocityField32    = uw.mesh.MeshVariable( mesh=mesh32,         nodeDofCount=dim )

    temperatureField32.load(inputPath+'Arrhenius_32_T.inp')
    velocityField32.load(inputPath+'Arrhenius_32_v.inp')
    pressureField32.load(inputPath+'Arrhenius_32_P.inp')
    
    temperatureField.data[:] = temperatureField32.evaluate(mesh)
    pressureField.data[:] = pressureField32.evaluate(mesh.subMesh)
    velocityField.data[:] = velocityField32.evaluate(mesh)
    # load summary statistics into arrays
    data = np.loadtxt(inputPath+'ArrSumary.txt', unpack=True )
    timeVal, vrmsVal, nuVal = data[0].tolist(), data[1].tolist(), data[2].tolist()<jupyter_output><empty_output><jupyter_text>**If _not_ loading from file: Initialise data**Start with a perturbed temperature gradient to speed up the convergence to the benchmark steady state solution.<jupyter_code>if(LoadFromFile == False):
    velocityField.data[:] = [0.,0.]
    pressureField.data[:] = 0.
    temperatureField.data[:] = 0.
    pertStrength = 0.1
    deltaTemp = tempMax - tempMin
    for index, coord in enumerate(mesh.data):
        pertCoeff = math.cos( math.pi * coord[0]/boxLength ) * math.sin( math.pi * coord[1]/boxLength )
        temperatureField.data[index] = tempMin + deltaTemp*(boxHeight - coord[1]) + pertStrength * pertCoeff
        temperatureField.data[index] = max(tempMin, min(tempMax, temperatureField.data[index]))
# initialise summary statistics arrays
    timeVal = []
    vrmsVal = []
    nuVal = []<jupyter_output><empty_output><jupyter_text>**Boundary conditions**This step is to ensure that the temperature boundary conditions are satisfied, as the initial conditions above may have been set to different values on the boundaries.<jupyter_code>for index in mesh.specialSets["MinJ_VertexSet"]:
    temperatureField.data[index] = tempMax
for index in mesh.specialSets["MaxJ_VertexSet"]:
    temperatureField.data[index] = tempMin<jupyter_output><empty_output><jupyter_text>**Conditions on the boundaries**Construct sets for the both horizontal and vertical walls. Combine the sets of vertices to make the I (left and right side walls) and J (top and bottom walls) sets. Note that both sets contain the corners of the box.<jupyter_code>iWalls = mesh.specialSets["MinI_VertexSet"] + mesh.specialSets["MaxI_VertexSet"]
jWalls = mesh.specialSets["MinJ_VertexSet"] + mesh.specialSets["MaxJ_VertexSet"]
allWalls = iWalls + jWalls

freeslipBC = uw.conditions.DirichletCondition( variable      = velocityField, 
                                               indexSetsPerDof = (allWalls,jWalls) )
tempBC     = uw.conditions.DirichletCondition( variable      = temperatureField, 
                                               indexSetsPerDof = (jWalls,) )<jupyter_output><empty_output><jupyter_text>**Plot initial conditions**Automatically scale the size of the vector arrows for the velocity field maximum.<jupyter_code>velmax = np.amax(velocityField.data[:])
if(velmax==0.0): velmax = 1.0
figtemp = glucifer.Figure()
tempminmax = fn.view.min_max(temperatureField)
figtemp.append( glucifer.objects.Surface(mesh, tempminmax) )
figtemp.append( glucifer.objects.VectorArrows(mesh, velocityField, scaling=0.1/velmax, arrowHead=0.2 ) )
figtemp.show()<jupyter_output><empty_output><jupyter_text>Set up material parameters and functions-----<jupyter_code># Rheology

delta_Eta = 1.0e6
activationEnergy = np.log(delta_Eta)
Ra = 1.e6
fn_viscosity = delta_Eta * fn.math.exp( - activationEnergy * temperatureField )

densityFn = Ra*temperatureField
# Define our vertical unit vector using a python tuple (this will be automatically converted to a function)
z_hat = ( 0.0, 1.0 )
# Now create a buoyancy force fn using the density.
buoyancyFn = z_hat*densityFn<jupyter_output><empty_output><jupyter_text>**Plot the viscosity**<jupyter_code>figEta = glucifer.Figure()
figEta.append( glucifer.objects.Surface(mesh, fn_viscosity) )
figEta.show()<jupyter_output><empty_output><jupyter_text>System setup-----**Setup a Stokes system: advanced solver settings**Set up parameters for the Stokes system solver. For PIC style integration, we include a swarm for the a PIC integration swarm is generated within. For gauss integration, simple do not include the swarm. Nearest neighbour is used where required.<jupyter_code>stokes = uw.systems.Stokes(velocityField=velocityField, 
                              pressureField=pressureField,
                              conditions=[freeslipBC,],
                              fn_viscosity=fn.exception.SafeMaths(fn_viscosity), 
                              fn_bodyforce=buoyancyFn )

solver=uw.systems.Solver(stokes)
solver.set_inner_method("mumps")
solver.set_penalty(1.0e7)<jupyter_output><empty_output><jupyter_text>**Create an advective diffusive system**<jupyter_code>advDiff = uw.systems.AdvectionDiffusion( temperatureField, tempDotField, velocityField, fn_diffusivity=1., conditions=[tempBC,], _allow_non_q1=True )<jupyter_output><empty_output><jupyter_text>Analysis tools-----**RMS velocity**Set up integrals used to calculate the RMS velocity.<jupyter_code>v2sum_integral  = uw.utils.Integral( mesh=mesh, fn=fn.math.dot(velocityField, velocityField) ) 
volume_integral = uw.utils.Integral( mesh=mesh, fn=1. )<jupyter_output><empty_output><jupyter_text>**Nusselt number**<jupyter_code>def FindNusseltNumber(temperatureField, linearMesh, xmax, zmax):
    tempgradField = temperatureField.fn_gradient
    vertGradField = tempgradField[1]
    BottomInt = 0.0
    GradValues = temperatureField.fn_gradient[1].evaluate(mesh.specialSets["MaxJ_VertexSet"])
    TopInt = sum(GradValues)
    for index in mesh.specialSets["MinJ_VertexSet"]:
        BottomInt += temperatureField.data[index]
    Nu = -zmax*TopInt/BottomInt
    return Nu[0]<jupyter_output><empty_output><jupyter_text>Main simulation loop-----The main time stepping loop begins here. Before this the time and timestep are initialised to zero and the output statistics arrays are set up. Since this may be a continuation of the saved data with associated summary statistics then check if there are existing time values first. If there is existing data then add new simulation statistics after the existing data.<jupyter_code>steps_prev = len(timeVal)
steps = 0
steps_end = 10000
step_out = 25

# Set time to zero, unless we are loading from file.
try:
    time_start = timeVal[-1]
except:
    time_start = 0.0
print 'Begining at t = ',time_start,' after having completed ',steps_prev,' steps'

simtime = time_start
# Setup clock to calculate simulation CPU time.
start = time.clock()

# Perform steps_end steps.
rms_v, nu_no = 0.0, 0.0
while steps<steps_end:
    # Get solution for initial configuration.
    solver.solve()
    # Retrieve the maximum possible timestep for the AD system.
    dt = advDiff.get_max_dt()
    if steps == 0:
        dt = 0.
    # Advect using this timestep size.
    advDiff.integrate(dt)
    
    # Calculate the RMS velocity.
    rms_v = math.sqrt(v2sum_integral.evaluate()[0]/volume_integral.evaluate()[0])
    nu_no = FindNusseltNumber(temperatureField, mesh, boxLength, boxHeight)
    
    # Increment time and store results.
    simtime += dt
    steps += 1
    vrmsVal.append( rms_v )
    timeVal.append( simtime )
    nuVal.append(nu_no)
                   
    if steps ==0 or steps % step_out == 0:
        print "steps = {:04d}; time = {:.6e}; vrms = {:6.2f}; nusselt = {:6.4f}; CPU = {:4.1f}s".format(
               steps, simtime, rms_v, nu_no, time.clock()-start)
        
v2sum = v2sum_integral.evaluate()
rms_v = rms_v = math.sqrt(v2sum_integral.evaluate()[0]/volume_integral.evaluate()[0])
# nu_no = -surface_Tgradient.integrate()[0] / basalT.integrate()[0]
nu_no = FindNusseltNumber(temperatureField, mesh, boxLength, boxHeight)

print "steps = {:04d}; time = {:.6e}; vrms = {:6.2f}; nusselt = {:6.4f}; CPU = {:4.1f}s".format(
               steps, simtime, rms_v, nu_no, time.clock()-start)<jupyter_output>steps = 0025; time = 2.929688e-03; vrms =   2.69; nusselt = 1.0000; CPU = 19.0s<jupyter_text>Post simulation analysis-----**Check CPU timing**Output timing for calculation per simulation time step<jupyter_code>if(steps!=steps_prev):
    avtime = (time.clock() - start)/float(steps)
else:
    avtime = 0.0

print "Average time per timestep = ",avtime," seconds over ", steps, " steps"<jupyter_output><empty_output><jupyter_text>Pre-run model took about $1.74$ seconds per timestep and was run for $13000$ timesteps for $Ra = 10^6$ for a $64\times64$ grid. **Save data to files**Save system summary data; $v_{rms}$ and the Nusselt number against time.<jupyter_code>np.savetxt( outputPath+'MandS_Summary.out', np.c_[timeVal, vrmsVal, nuVal], header="Time, VRMS, Nusselt" )<jupyter_output><empty_output><jupyter_text>Save final temperature, velocity and pressure fields.<jupyter_code>temperatureField.save(outputPath+'MandS_T.out')
velocityField.save(outputPath+'MandS_v.out')
pressureField.save(outputPath+'MandS_p.out')<jupyter_output><empty_output><jupyter_text>**Plot temperature field with velocity vectors overlaid**<jupyter_code>velmax = np.amax(velocityField.data[:])
figVT = glucifer.Figure()
figVT.append( glucifer.objects.Surface(mesh, temperatureField) )
figVT.append( glucifer.objects.VectorArrows(mesh, velocityField, scaling=0.1/velmax, arrowHead=0.2) )
figVT.show()<jupyter_output><empty_output><jupyter_text>**Plot system summary information**If the initial conditions were loaded up from a data file then these arrays will contain the saved data as well as results from the simulation just run.<jupyter_code>pylab.rcParams[ 'figure.figsize'] = 14, 8
pyplot.figure(1)
pyplot.subplot(211)
pyplot.plot(timeVal, nuVal)
pyplot.ylabel('Nusselt Number')

pyplot.subplot(212)
pyplot.plot(timeVal,vrmsVal)
pyplot.xlabel('Time')
pyplot.ylabel('Vrms')
pyplot.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10860271506787669, 'fraction_numerical': 0.019275481887047176, 'mean_word_length': 4.31870761866773, 'pattern_counts': {'":': 0, '<': 94, '<?xml version=': 0, '>': 93, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '22', 'dir': 'jupyter-structured-clean-dedup', 'id': '514022', 'license': 'MIT', 'n_tokens_mistral': 4366, 'n_tokens_neox': 4129, 'n_words': 1085, 'path': 'notebooks/publications/MoresiSolomatov-1995/MoresiAndSolomatov1995.ipynb', 'repo_name': 'lmoresi/uw-docker-landing-page', 'size': '13369'}