Hi all,
I tried to calculate a subduction model using the particle method, which is my initial setup
###### Initial compositional fields #####
subsection Compositional fields
** set Number of fields = 12**
** set Names of fields =Mantle, ContinentalUpperCrust, ContinentalLowerCrust,ContinentalMantle, CratonMantle , OceanicCrust,OceanicMantle, sediment, sediment_age , noninitial_plastic_strain, plastic_strain , viscous_strain**
** set Compositional field methods = particles, particles,particles, particles,particles, particles,particles, particles,particles, particles,particles, particles**
** set Mapped particle properties = Mantle:initial Mantle , ContinentalUpperCrust:initial ContinentalUpperCrust , ContinentalLowerCrust:initial ContinentalLowerCrust , ContinentalMantle:initial ContinentalMantle , CratonMantle:initial CratonMantle , OceanicCrust:initial OceanicCrust , OceanicMantle:initial OceanicMantle , sediment:initial sediment , sediment_age:initial sediment_age , noninitial_plastic_strain:initial noninitial_plastic_strain , plastic_strain:initial plastic_strain , viscous_strain:initial viscous_strain**
end
subsection Initial composition model
** set List of model names = ascii data , fault noise**
** set List of model operators = add**
** subsection Ascii data model**
** set Data directory = /work/home/hanx23/aspect/aspect_precipitation/aspect/workspace/test/**
** set Data file name = Composition.txt**
** end**
** subsection Fault Noise**
** set Standard deviation of Gaussian noise amplitude distribution = 200e3 **
** set Maximum amplitude of Gaussian noise amplitude distribution = 0.5**
** set Random number generator seed = 631**
** set Depth around which Gaussian noise is smoothed out = 40e3 **
** set Halfwidth with which Gaussian noise is smoothed out in depth = 10e3 **
** # The rift polygon specified in cartesian surface coordinates (x and y)**
** set Fault axis line segments = 1350e3**
** # The resolution of the table that holds the noise values**
** # Bilinear interpolation between these values is performed when necessary**
** set Grid intervals for noise X or radius = 2880**
** set Grid intervals for noise Y or longitude = 1280**
** end**
end
###### Initial temperature model ######
subsection Initial temperature model
** set List of model names = ascii data , adiabatic**
** set List of model operators = add, minimum**
** subsection Ascii data model**
** set Data directory = /work/home/hanx23/aspect/aspect_precipitation/aspect/workspace/test/**
** set Data file name = Temperature.txt**
** end**
** subsection Adiabatic**
** subsection Function**
** set Variable names = t**
** set Function expression = if(t>=4000,1,0);**
** if(t>=4000,1,0);**
** if(t>=4000,1,0);**
** if(t>=4000,1,0);**
** if(t>=4000,1,0);**
** if(t>=4000,1,0);**
** if(t>=4000,1,0);**
** if(t>=4000,1,0);**
** if(t>=4000,1,0);**
** if(t>=4000,1,0);**
** if(t>=4000,1,0);**
** if(t>=4000,1,0)**
** end**
** set Age top boundary layer = 0**
** set Age bottom boundary layer = 0**
** end**
end
####### Boundary conditions #######
####### Boundary composition model #######
subsection Boundary composition model
** set List of model names = initial composition**
** set Fixed composition boundary indicators = right ,bottom**
** #set Allow fixed composition on outflow boundaries = true**
** subsection Initial composition**
** end**
end
####### Boundary temperature model #######
subsection Boundary temperature model
** set Fixed temperature boundary indicators = top , bottom , left , right**
** set List of model names = initial temperature**
end
####### Boundary velocity model ######
####### Right boundary velocity 2 cm/yr Bottom velocity #######
subsection Boundary velocity model
** set Prescribed velocity boundary indicators = right x:function , bottom y:function**
** subsection Function**
** set Variable names = x,y**
** set Function constants = cm=0.01 , yr = 1 , VelContinental = 2 , Ay = 900e3 , HFore = 130e3 , XLong = 2000e3**
** set Function expression =if( y>=(Ay - HFore) , -VelContinentalcm/yr ,*
** if( y<(Ay - HFore) && y>=(Ay - HFore - HFore) , -VelContinental/(HFore)*(y - (Ay - HFore - HFore ) )cm/yr , 0 ) );\ **
** 0*
** end**
end
####### Boundary traction model #######
####### All the left and bottom are open #######
subsection Boundary traction model
** set Prescribed traction boundary indicators = left x : initial lithostatic pressure **
** subsection Initial lithostatic pressure**
** set Number of integration points = 1800**
** set Representative point = 150000, 900000**
** end**
end
######## Coupling Fastscape ########
subsection Mesh deformation
** set Mesh deformation boundary indicators = top : fastscape**
** set Additional tangential mesh velocity boundary indicators = left, right**
** # Renamed free surface stabilization and moved it into general mesh deformation section.**
** set Surface stabilization theta = 0.5**
** subsection Fastscape**
** set Vertical exaggeration = 1**
** set Maximum timestep = 5000**
** set Number of steps = 10**
** set Surface resolution = 10**
** set Resolution difference = 4**
** set Additional fastscape refinement = 0**
** set Use marine component = true**
** set Use ghost nodes = true**
** set Y extent in 2d = 50e3**
** set Sediment rain = 0,0 # m/yr**
** set Sediment rain intervals = 0**
** set Use velocities = true**
** set Use precipitation model = true**
** subsection Boundary conditions**
** set Bottom = 1**
** set Right = 1**
** set Top = 1**
** set Left = 1**
** end**
** subsection Erosional parameters**
** set Drainage area exponent = 0.4 #m**
** set Slope exponent = 1 #n**
** set Multi-direction slope exponent = -1 #p**
** set Bedrock diffusivity = 5e-3 #kd**
** set Bedrock river incision rate = 1e-5 #kf**
** set Bedrock deposition coefficient = 1 #G**
** set Sediment diffusivity = -1**
** set Sediment river incision rate = -1 **
** set Sediment deposition coefficient = 1 #G**
** end**
** subsection Marine parameters**
** set Sea level = 3000**
** set Sand porosity = 0**
** set Shale porosity = 0**
** set Sand e-folding depth = 0**
** set Shale e-folding depth = 0**
** set Sand-shale ratio = 1**
** set Depth averaging thickness = 1e2**
** set Sand transport coefficient = 200**
** set Shale transport coefficient = 200**
** end**
** end**
** end**
####### Postprocess #######
subsection Postprocess
set List of postprocessors = velocity statistics, basic statistics, temperature statistics, heat flux statistics, visualization, mass flux statistics, composition statistics, topography, Stokes residual, particles
subsection Visualization
set Interpolate output = true
set Write higher order output = true
set List of output variables = viscosity, heat flux map, vertical heat flux, density, strain rate, depth, named additional outputs, maximum horizontal compressive stress, principal stress
set Time between graphical output = 0
set Number of grouped files = 1
set Point-wise stress and strain = true
end
subsection Particles
** set Number of particles = 1000000**
** set Load balancing strategy = remove and add particles**
** set Time between data output = 0**
** set Minimum particles per cell = 10**
** set Maximum particles per cell = 50**
** set Data output format = vtu**
** set List of particle properties = initial composition, velocity**
** set Particle generator name = random uniform**
** set Update ghost particles = true**
** end**
end
This is my output file:
output.txt (47.6 KB)
I think that excessive grid deformation is causing the issue. To address this, I reduced the particle generation boundaries to make it smaller than the grid and increased the number of iterations for particle generation, aiming to improve the chances that the particles will fall within the grid.
But the problem remains：
output2.txt (47.5 KB)
I try to ignore this error and let the program continue, I want to see where is the mesh deformation that prevents it from generating particles.
I see that the right boundary grid is getting rid of particles.
I set a velocity boundary condition of 2 cm/year on the right boundary, where the velocity enters the domain. The mesh doesn’t appear to be significantly deformed, so I don’t understand why particles cannot be generated. I would appreciate any help or insight.
Best
Xiao
