[25806] | 1 | #Test Name: TransientFrictionSchoof
|
---|
| 2 | import numpy as np
|
---|
| 3 |
|
---|
| 4 | from frictionschoof import frictionschoof
|
---|
[26638] | 5 | from socket import gethostname
|
---|
[25806] | 6 | from model import *
|
---|
| 7 | from parameterize import parameterize
|
---|
| 8 | from setflowequation import setflowequation
|
---|
| 9 | from setmask import setmask
|
---|
| 10 | from solve import solve
|
---|
| 11 | from transient import transient
|
---|
| 12 | from triangle import triangle
|
---|
| 13 |
|
---|
| 14 |
|
---|
| 15 | md = triangle(model(), '../Exp/Square.exp', 200000.)
|
---|
| 16 | md = setmask(md, '', '')
|
---|
| 17 | md = parameterize(md, '../Par/SquareSheetConstrained.py')
|
---|
| 18 | md = setflowequation(md, 'SSA', 'all')
|
---|
| 19 |
|
---|
| 20 | # Use Schoof's law
|
---|
| 21 | Cmax = 0.8
|
---|
| 22 | md.friction = frictionschoof()
|
---|
| 23 | md.friction.m = 1.0 / 3.0 * np.ones((md.mesh.numberofelements, 1))
|
---|
| 24 | md.friction.Cmax = Cmax * np.ones((md.mesh.numberofvertices, 1))
|
---|
| 25 | md.friction.C = 200 * np.ones((md.mesh.numberofvertices, 1))
|
---|
| 26 |
|
---|
| 27 | # Control parameters
|
---|
| 28 | md.inversion.iscontrol = 1
|
---|
| 29 | md.inversion.control_parameters = ['FrictionC']
|
---|
| 30 | md.inversion.min_parameters = 1. * np.ones((md.mesh.numberofvertices, 1))
|
---|
| 31 | md.inversion.max_parameters = 10000. * np.ones((md.mesh.numberofvertices, 1))
|
---|
| 32 | md.inversion.nsteps = 2
|
---|
| 33 | md.inversion.cost_functions = [102, 501]
|
---|
| 34 | md.inversion.cost_functions_coefficients = np.ones((md.mesh.numberofvertices, 2))
|
---|
| 35 | md.inversion.cost_functions_coefficients[:, 1] = 2e-7
|
---|
| 36 | md.inversion.gradient_scaling = 3. * np.ones((md.inversion.nsteps, 1))
|
---|
| 37 | md.inversion.maxiter_per_step = 2 * np.ones((md.inversion.nsteps, 1))
|
---|
| 38 | md.inversion.step_threshold = 0.3 * np.ones((md.inversion.nsteps, 1))
|
---|
| 39 | md.inversion.vx_obs = md.initialization.vx
|
---|
| 40 | md.inversion.vy_obs= md.initialization.vy
|
---|
| 41 |
|
---|
[26638] | 42 | md.cluster = generic('name', gethostname(), 'np', 3)
|
---|
[25806] | 43 | md = solve(md, 'Stressbalance')
|
---|
| 44 |
|
---|
| 45 | #Fields and tolerances to track changes
|
---|
| 46 | field_names = ['Gradient', 'Misfits', 'FrictionC', 'Pressure', 'Vel', 'Vx', 'Vy']
|
---|
| 47 | field_tolerances = [1e-12, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13]
|
---|
| 48 | field_values = [
|
---|
| 49 | md.results.StressbalanceSolution.Gradient1,
|
---|
| 50 | md.results.StressbalanceSolution.J,
|
---|
| 51 | md.results.StressbalanceSolution.FrictionC,
|
---|
| 52 | md.results.StressbalanceSolution.Pressure,
|
---|
| 53 | md.results.StressbalanceSolution.Vel,
|
---|
| 54 | md.results.StressbalanceSolution.Vx,
|
---|
| 55 | md.results.StressbalanceSolution.Vy
|
---|
| 56 | ]
|
---|