I tried to generate plate-driven flow through the 3-D chunk geometry.
I used 'Boundary velocity model" with Prescribed velocity boundary indicators (outer) and Tangential velocity boundary indicator (inner), like below:
subsection Boundary velocity model
set Tangential velocity boundary indicators = inner
set Prescribed velocity boundary indicators = outer:function
subsection Function
set Coordinate system = spherical
set Use spherical unit vectors = true
set Variable names = r, phi, theta
set Function constants = mm=0.001, year=1, km=1000, nver=-14.02, evel=21.27
set Function expression = 0; if(r>=6330km, evelmm/year, 0); if(r>=6330km, nvelmm/year, 0)
end
end
In the description of “Use spherical unit vectors”, Positive values of first, second, and third components indicates up, east, and north.
I expected that the NE direction flow was generated, because east component was positive and north component was negative.
However, SE direction flow was generated. Could you please any comments about this isssue?
My initial temperature model consists of spherical coordinates (r, phi, theta).
I assigned that r means earth_radius-depth, phi means longitude, theta means co-latitude (i.e., 90 - latitude).
If I was correct, I think positive third component of “Use spherical unit vectors” should indicate south direction.
The radius, longitude, co-latitude coordinates for ascii data of initial temperature model is correct?
Can you send over the entire prm file and a picture of the flow field?
If possible, can you make the PRM file low resolution so it is easy for us to re-run the model on one or two cores?
What are the boundary conditions on the side walls? By default, these boundaries will be assigned pressure boundary conditions if nothing is specified in the PRM file, and this may dominate the flow field. I would first make sure you have the desired boundary condition (velocity or lithostatic pressure) on these boundaries.
I have attached my entire prm file here. original.prm (3.4 KB)
I would like to provide low resolition prm file, I’ll think about how to modify the prm file and provide it to you.
I previously assigned default boundaries because when I assigned free slip to the side walls, simulation did not finish.
When free slip boundaries are applied to the side walls, the simulation does not proceeded to the next step in the image below (I waited it more than one week, it did not proceed next step).
@kjun825 - Given you are waiting some time for the current models to converge, I think the most practical solution is to develop a much smaller (lower resolution) model that can be run quickly for testing.
From the other forum posts regarding the prescribed velocities, it sounds like that may also lead to improved convergence behavior.
Let us know if we can continue to help debug after designing a lower resolution test case.
This prm file runs quickly on a single core due to the low viscosity and large grid size. low_visco_large_grid.prm (3.8 KB)
In the high resolution simulatation with free slip boundary conditions, it failed to converge (solving stokes system did not finished). high_resolution_with_freeslip.prm (3.5 KB)
However, when I used an linitial lithostatic pressure boundary condition in the high-resolution case, it did converge. high_resolution_with_traction.prm (3.8 KB)
I aimed to generate NW to SE direction of surface motion, but the simulation produced surface motion in the SW to NE direction of surface motion.
In the high resolution simulatation with free slip boundary conditions, it failed to converge (solving stokes system did not finished).
My expectation would actually be that this model should not converge.
If you are applying a prescribed velocity (plate motion) at the surface, this inherently produces a velocity discontinuity when the surface intersects the side walls.
Using a lithostatic pressure condition on the side walls allows material beneath the surface to flow across the boundaries in response to the prescribed surface velocity and any additional driving forces in your model domain.
Can you provide an overview of what your are trying to model with this setup? Knowing the end goal would allows to provide more helpful suggestions.
My goal is to simulate the interaction between plate-driven mantle flow and lithosphere of varing thickness.
To observe this interaction, I attemted to construct variations in lithosphric thickness at the top and generat the NW-SE direction of flow (but it generaged SW-NE direction of flow).
The reason I assigned a free-slip boundary condition is that I want to focus on the interaction around the lithosphere without effects from outside the domain.
If setting a lithostatic preesure condition on the side wall is more appropriate, I can ust that as the boundary condition.
I have a one more question… In the forum people usually set the bottom-right point as a representative point for initial lithostatic pressure.
In my case, should be it (5710000, 150, 13) or (6370000, 150, 13)? I assume the radius of the Earth is 6370000.
If setting a lithostatic preesure condition on the side wall is more appropriate, I can just that as the boundary condition.
This is a really a model design question that you must ultimately make a decision on, but here are what I view as feasible options for the side wall boundary conditions
Use lithostatic pressure boundary condition
Apply the surface velocity function to the side walls, but have the velocities gradate in some fashion from the surface to the base of the model where a no slip condition is applied.
Use a chunk geometry with lithosphere boundary indicators to have separate boundary conditions in the lithosphere and deeper parts of the model domain.
I have a one more question… In the forum people usually set the bottom-right point as a representative point for initial lithostatic pressure.
In my case, should be it (5710000, 150, 13) or (6370000, 150, 13)? I assume the radius of the Earth is 6370000.
Per the documentation here, the depth/radius value actually does not matter. You should need to pick the desired lat/lon value for a spherical geometry.
I’ll try construct the model again, refering to your comments.
Lastly, I assigned a positive value to the east component and a negative value to the north component in “Prescribed velocity boundary indicators” section to generate the SE direction of the surface flow.
According to the manual, if I set “Use spherical unit vector = true”, then positive values indicate the east and north directions.
However, when I assign a negative value to the north component (the south direction), the flow moves toward the north (the NE direction).
Could you please provide any comment on this? I don’t rule out the possibility that I may have misunderstood the manual.
We changed the direction of the “N-S vector” 5 years ago (new vs old) and must have forgotten to change the manual.
For others trying to figure out what the manual should say, the three numbers passed correspond to the velocity components in the local “Up”, “East” and “South” directions (in that order).
Thanks for pointing this out to us, we’ll fix the manual (you can track the issue here).
