#!/usr/bin/env python # coding: utf-8 ### Arc-continent collision Bdls - script - 10-03-22 import UWGeodynamics as GEO import numpy as np from underworld import function as fn from UWGeodynamics import visualisation as glucifer from MechanicalProperties import PlateProperties from ModelGeometry import SubductionCreator,interpolateTracer,rmRepeated,fuseListM,ListToNd,generateWeakzone import os import underworld as uw # %% # shortcuts for parallel wrappers barrier = GEO.uw.mpi.barrier rank = GEO.rank #Units u = GEO.UnitRegistry ######################################## #Rebecca et al., Scaling dRho = 80. * u.kilogram / u.meter**3 # matprop.ref_density g = 9.8 * u.meter / u.second**2 # modprop.gravity #H = 800. * u.kilometer # modprop.boxHeight-- OR 840 dependig if I include the sticky air layer H= 900. * u.kilometer # lithostatic pressure for mass-time-length ref_stress = dRho * g * H # viscosity of upper mante for mass-time-length ref_viscosity = 1e20 * u.pascal * u.seconds #References ref_time = ref_viscosity/ref_stress ref_length = H ref_mass = (ref_viscosity*ref_length*ref_time) #ref_temperature = modprop.Tint - modprop.Tsurf KL = ref_length KM = ref_mass Kt = ref_time #KT = ref_temperature GEO.scaling_coefficients["[length]"] = KL GEO.scaling_coefficients["[time]"] = Kt GEO.scaling_coefficients["[mass]"]= KM #GEO.scaling_coefficients["[temperature]"] = KT ########################################## ######################################################### ## These are the C parameters of Crameri et al., 2011 for choosing a proper sticky air layer thickness and viscosity H_sticky=200 #km # %% # #Sticky air calcualtion - Stokes L=900#980 c=(L/H_sticky)**3*(1e20/1e20)*(1/16)*(80/3400) c # %% # #Sticky air calcualtion - Isostacy L=2600#980 c=(L/H_sticky)**3*(1e20/1e20)*(3/(16*(np.pi**3))) c ##################################################### #Model creation res=(650,300) ## Resolution. This gives 5 km in x and 3 km in Y #res=(300,250) #Model Dimensions Model = GEO.Model(elementRes=res, #296,160 #296,128, #120,64 minCoord=(-1300.0 * u.kilometer, -700.0 * u.kilometer), maxCoord=(1300.0* u.kilometer, 200.0 * u.kilometer), gravity=(0.0, -9.81 * u.meter / u.second**2)) #Output folder definition outputPath="Arc_collision_bdls_2e-5" Model.outputDir =outputPath # %% outputPath = os.path.join(os.path.abspath("."), outputPath) if rank==0: if not os.path.exists(outputPath): os.makedirs(outputPath) barrier() #Model Geometry #Decoupling Layer- 70 Km for 296, for 120, 110 DipAngle=30 DipLen=300 #original SP length=2400 #SubductionCreator(Model,y0,thickness, dipAngle, dipLength, maxLength, orientation, SLayers,OLayers, ExLens,bStrips,tD=False) geometry=SubductionCreator(Model,0,100,DipAngle,DipLen,1700,-1,4,2,(210,210,210,210),(40,40,40,40),(True,True, True, True),0) xcoords = GEO.uw.function.input()[0] ycoords = GEO.uw.function.input()[1] #Model limits orientation=geometry[6][0] xlimit=GEO.nd(geometry[6][1]) #cratLim=GEO.nd(geometry[8][2]) ## Calculating material properties for all materials in the model. #PlateProperties(Nlayers,crustThickness,crustDensity,mantleDensity,plateThickness,oceanic,age, cohesion,friction,arc,arcDensity,arcThickness,depthToMantle) #Oceanic Plate data1=PlateProperties(4,7.,2900.,3400.,100.,True,80.,12.5,12.5,0.066,0.033,False,0.,0.,100.,2.) #Cratonic continental - Values can be +5MPa!! data2=PlateProperties(2,40.,2700.,3400.,150.,False,80.,15.,10.,0.15,0.08,False,0.,0.,150.,.2) #Back-arc data3=PlateProperties(2,20.,2800.,3400.,100.,False,80.,10.,7.5,0.005,0.0035,False,0.,0.,80.,2.) #the last specifies how much was extended the lithosphere #Arc Crust data4=PlateProperties(4,7.,2900.,3400.,100.,True,80.,12.5,12.5,0.066,0.033,True,2838.,25.,100.,2.) #eclogite-Properties data5=PlateProperties(4,2.,3450.,3400.,100.,True,80.,12.5,12.5,0.066,0.033,False,0.,0.,100.,2.) #Transitional crust data6=PlateProperties(2,25.,2700.,3400.,100.,False,80.,12.5,12.5,0.066,0.033,False,0.,0.,100.,2.) # %% #Densities without units refD=3400.* u.kilogram / u.metre**3 #Mantle Density #(Subducting plate) l1d=data1[2][0]* u.kilogram / u.metre**3 l2d=data1[2][1]* u.kilogram / u.metre**3 l3d=data1[2][2]* u.kilogram / u.metre**3 l4d=data1[2][3]* u.kilogram / u.metre**3 #Cratonic overriding plate cl1=data2[2][0]* u.kilogram / u.metre**3 cl2=data2[2][1]* u.kilogram / u.metre**3 #Back-Arc overriding plate bl1=data3[2][0]* u.kilogram / u.metre**3 bl2=data3[2][1]* u.kilogram / u.metre**3 #Arc density al1=data4[2][0]* u.kilogram / u.metre**3 al2=data4[2][1]* u.kilogram / u.metre**3 #Oceanic crust to Eclogite e1d=data5[2][0]* u.kilogram / u.metre**3 #e1d= 3450* u.kilogram / u.metre**3 #e1d=3363.51575* u.kilogram / u.metre**3 #Lower Mantle density - 1.6 +-1 to 5+-2 percentage the upper mantle density - (1.6+7)/2/100 - 3400.+(3400*0.043) #LMdensity=3546.2* u.kilogram / u.metre**3 LMdensity=3488.4* u.kilogram / u.metre**3 #transition_crust tra1=data6[2][0]* u.kilogram / u.metre**3 tra2=data6[2][1]* u.kilogram / u.metre**3 ## Generate a weak zone that will enable subduction initiation weakShape=generateWeakzone(geometry[6][1],0*u.kilometer, 30,DipAngle, DipLen) #This cointains an initial plastic damage of 0.2 weakShape2=generateWeakzone(geometry[6][1],0*u.kilometer, 250,DipAngle, DipLen) #This allows plastic strain around the trench ## Calculate Subduction geometry Geometry #SubductionCreator(Model,y0,thickness, dipAngle, dipLength, maxLength, orientation, SLayers,OLayers, ExLens,bStrips,tD=False) #return subducting,overriding, weak,arc,bstop #Add materials to the model Air_shape=GEO.shapes.Layer(top=200.*u.kilometer, bottom=0.*u.kilometer) Air=Model.add_material(name="Air", shape=Air_shape) UMantle =Model.add_material(name="UpperMantle", shape=GEO.shapes.Layer(top=0.*u.kilometer, bottom=-700.*u.kilometer)) OLithosphere1F = Model.add_material(name="SubductingPlateL1", shape=geometry[0][2][0]) OLithosphere2F=Model.add_material(name="SubductingPlateL2", shape=geometry[0][2][1]) OLithosphere3F=Model.add_material(name="SubductingPlateL3", shape=geometry[0][2][2]) OLithosphere4F=Model.add_material(name="SubductingPlateL4", shape=geometry[0][2][3]) Clithosphere1= Model.add_material(name="Overriding plate Crust", shape=geometry[1][0]) Clithosphere2= Model.add_material(name="Overriding plate Lithosphere", shape=geometry[1][1]) Clithosphereweak1=Model.add_material(name="Overriding plate Weak Crust", shape=geometry[2][0]) Clithosphereweak2=Model.add_material(name="Overriding plate Weak Lithosphere", shape=geometry[2][1]) OArc1=Model.add_material(name="IntraOceanicArc", shape=geometry[3][0]) Ctransitional1=Model.add_material(name="OP_transitional1", shape=geometry[12][0]) Ctransitional2=Model.add_material(name="OP_transitional2", shape=geometry[12][1]) sediment = Model.add_material(name="Sediment") # The shape objetcs below will be used to remove plastic strain from model walls. This causes plastic strain to localize in the collisional zone, where it should #Strain-backstop OP plate shapes back1=geometry[11][0] back2=geometry[11][1] back3=geometry[11][2] back4=geometry[11][3] #Transitional crust shapes ctra1=geometry[12][0] ctra2=geometry[12][1] #cratonic continental plate shapes cratS1=geometry[1][0] cratS2=geometry[1][1] # %% #Preview of 2D materials-Materials Field (from swarm) - Glucifer plot Fig = glucifer.Figure(figsize=(1200,400)) Fig.Points(Model.swarm, Model.materialField,fn_size=2.0, discrete=True) #Fig.Surface(Model.mesh,Model.projMaterialField,fn_size=2.0) Fig.show() # %% #eclogite density fn conditions=[(Model.y > GEO.nd(-150.*u.kilometer),GEO.nd(l1d)), (True,GEO.nd(e1d)), ] eclogiteFn=fn.branching.conditional(conditions) #avoid continental crust reaching --400 km to come back surface conditionsSubCrust=[(Model.y > GEO.nd(-400.*u.kilometer),GEO.nd(bl1)), (True,GEO.nd(refD)), ] #avoid air-mantle contact - All mantle material reaching the surface becomes basalt conditionsAvoid=[(((Model.x > GEO.nd(0.*u.kilometer)) & (Model.y > GEO.nd(-5.*u.kilometer))),GEO.nd(l1d)), (True,GEO.nd(refD)), ] AvoidFn=fn.branching.conditional(conditionsAvoid) #all continental crust below 400 becmes mantle material to avoid chucks of continental crust to ascend through the mantle CrustFn=fn.branching.conditional(conditionsSubCrust) #initial plastic Strain - This function "damages" the SP-OP interface to facilitate subduction initiation plasticCondition=[(Model.y (dte / 3.): self._dt = dte / 3. comm.Barrier() _update(Model) #self._update() self.step += 1 self.stepDone += 1 self._ndtime += self._dt checkpointer.checkpoint() if rank == 0: string = """Step: {0:5d} Model Time: {1:6.1f} dt: {2:6.1f} ({3})\n""".format( self.stepDone, _adjust_time_units(self.time), _adjust_time_units(self._dt), datetime.now().strftime('%Y-%m-%d %H:%M:%S')) sys.stdout.write(string) sys.stdout.flush() self._post_solve() return 1 #####From UW Source Code. Beucher et al., 2022- this allows to run the model using voronoi swarm def stokes_SLE(Model): """ Stokes SLE """ self=Model if any([material.viscosity for material in self.materials]): conditions = list() conditions.append(self.velocityBCs) if self._stressBCs: conditions.append(self.stressBCs) fssa = None if self._fssa_factor: fssa = self._fssa self._stokes_SLE = uw.systems.Stokes( velocityField=self.velocityField, pressureField=self.pressureField, voronoi_swarm = self.swarm, conditions=conditions, _fn_fssa = fssa, fn_viscosity=self._viscosityFn, fn_bodyforce=self._buoyancyFn, fn_stresshistory=self._elastic_stressFn, fn_one_on_lambda=self._lambdaFn) return self._stokes_SLE Model._solver = uw.systems.Solver(stokes_SLE(Model)) #Solver Parameters Model.solver.set_inner_method("mumps") Model.solver.set_penalty(1e6) #GEO.rcParams["initial.nonlinear.tolerance"] = 1e-4 # %% #Fields to Save outputss=['temperature', 'pressureField', 'strainRateField', 'velocityField', 'projStressField', 'projTimeField', 'projMaterialField', 'projViscosityField', 'projStressField', 'projMeltField', 'projPlasticStrain', 'projDensityField', 'projStressTensor', ] GEO.rcParams['default.outputs']=outputss #This is needed just in case the model do not converges with the default configuration of the solver #solver = Model.solver #scr_rtol = 1e-6 ## Schur complement solver options #solver.options.scr.ksp_rtol = scr_rtol #solver.options.scr.ksp_type = "fgmres" ##solver.options.main.list() # ## Inner solve (velocity), A11 options #solver.options.A11.ksp_rtol = 1e-1 * scr_rtol #solver.options.A11.ksp_type = "fgmres" #Enables the coupling with badlands Model.surfaceProcesses = GEO.surfaceProcesses.Badlands(airIndex=[Air.index],sedimentIndex=sediment.index, XML="./ressources/badlands3.xml", resolution=4. * u.kilometer,restartFolder="outbdls", checkpoint_interval=0.2 * u.megayears,restartStep=165) def post_hook(): vrms = Model.stokes_SLE.velocity_rms() step = Model.step time = Model.time.m_as(u.year) if rank == 0: with open(fout,'a') as f: f.write(f"{step}\t{time:5e}\t{vrms:5e}\n") #This removes plastic damage from model left and right walls fact1=InfplasticFn.evaluate(Model.swarm) Model.plasticStrain.data[:] =Model.plasticStrain.data[:] * fact1 Model.post_solve_functions["B"] = post_hook #### This handles the initial instabilities in the model. We aim to "equilibrate the model" by removing all plastic strain during the first 1-3 Myr. ### We only allow plastic strain to localize in the trench during this first 1-3 Myr. ## After that the Equilibrium variable becomes True and the time step increases. MaUnit=1000000#*u.years #Equilibrium=False #Equilibrium=True restart=True if restart==False: #Small time step at the beginning - Because isostatic compesation occurs here run_for_adaptative(Model,duration=(0.2*MaUnit)*u.year,checkpoint_interval=(0.2*MaUnit)*u.year,dt=500*u.years,restartDir=outputPath) # All plastic damage duirng equilibration is removed Model.plasticStrain.data[:] = 0. #It gives some time the model to take "impulse" again as plastic damage was removed everywhere in the model run_for_adaptative(Model,duration=(0.2*MaUnit)*u.year,checkpoint_interval=(0.2*MaUnit)*u.year,dt=1500*u.years,restartDir=outputPath) #A little bigger time step during subduction initiation run_for_adaptative(Model,duration=(6.*MaUnit)*u.year,checkpoint_interval=(0.2*MaUnit)*u.year,dt=6000*u.years,restartDir=outputPath) #The maximum stable time step allowed by FSSA. run_for_adaptative(Model,duration=(25.*MaUnit)*u.year,checkpoint_interval=(0.2*MaUnit)*u.year,dt=7500*u.years,restartDir=outputPath) #Close to break-up the time step must decrease run_for_adaptative(Model,duration=(15.*MaUnit)*u.year,checkpoint_interval=(0.2*MaUnit)*u.year,dt=5000*u.years,restartDir=outputPath) else: run_for_adaptative(Model,duration=(20.*MaUnit)*u.year,checkpoint_interval=(0.2*MaUnit)*u.year,dt=5000*u.years,restartDir=outputPath,restartStep=164) run_for_adaptative(Model,duration=(10.*MaUnit)*u.year,checkpoint_interval=(0.2*MaUnit)*u.year,dt=5000*u.years,restartDir=outputPath)