Changeset 17469
- Timestamp:
- 03/18/14 16:32:17 (11 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/test/Par/ValleyGlacierShelf.par
r17467 r17469 1 1 %Start defining model parameters here 2 x = md.mesh.x 3 y = md.mesh.y 4 xmin = min(x) 5 xmax = max(x) 6 ymin = min(y) 7 ymax = max(y) 8 Lx = (xmax-xmin) 9 Ly = (ymax-ymin) 10 xm = (xmin+xmax)/2. 11 ym = (ymin+ymax)/2. 2 x = md.mesh.x; 3 y = md.mesh.y; 4 xmin = min(x); 5 xmax = max(x); 6 ymin = min(y); 7 ymax = max(y); 8 Lx = (xmax-xmin); 9 Ly = (ymax-ymin); 10 xm = (xmin+xmax)/2.; 11 ym = (ymin+ymax)/2.; 12 12 13 13 %Geometry: U-shaped valley in y direction 14 thk_center = 1000. 15 thk_margin = 0.5*thk_center 16 bmax = 0. 17 bmin = -thk_center*md.materials.rho_ice/md.materials.rho_water 14 thk_center = 1000.; 15 thk_margin = 0.5*thk_center; 16 bmax = 0.; 17 bmin = -thk_center*md.materials.rho_ice/md.materials.rho_water; 18 18 19 alpha = 2./3. 20 slope = 0.9*(bmin-bmax)*(x-xmin)/(Lx*alpha) + 0.1*(bmin-bmax)*(y-ymin)/(Ly) + bmax 21 md.geometry.surface = (thk_center+bmax) + slope 22 md.geometry.bed = bmax + slope + 4./Ly^2*(thk_center-thk_margin)*(y-ym).^2) 23 md.geometry.thickness = md.geometry.surface - md.geometry.bed 24 md.geometry.bathymetry = md.geometry.bed 19 alpha = 2./3.; 20 slope = 0.9*(bmin-bmax)*(x-xmin)/(Lx*alpha) + 0.1*(bmin-bmax)*(y-ymin)/(Ly) + bmax; 21 md.geometry.surface = (thk_center+bmax) + slope ; 22 md.geometry.bed = bmax + slope + 4./Ly^2*(thk_center-thk_margin)*(y-ym).^2); 23 md.geometry.thickness = md.geometry.surface - md.geometry.bed; 24 md.geometry.bathymetry = md.geometry.bed; 25 25 26 26 %Mask 27 md.mask.ice_levelset = x - alpha*Lx 28 md.mask.groundedice_levelset = ones(md.mesh.numberofvertices,1) 27 md.mask.ice_levelset = x - alpha*Lx; 28 md.mask.groundedice_levelset = ones(md.mesh.numberofvertices,1); 29 29 30 30 %Initial velocity 31 md.initialization.vx = zeros(md.mesh.numberofvertices,1) 32 md.initialization.vy = zeros(md.mesh.numberofvertices,1) 33 md.initialization.vz = zeros(md.mesh.numberofvertices,1) 34 md.initialization.pressure = zeros(md.mesh.numberofvertices,1) 31 md.initialization.vx = zeros(md.mesh.numberofvertices,1); 32 md.initialization.vy = zeros(md.mesh.numberofvertices,1); 33 md.initialization.vz = zeros(md.mesh.numberofvertices,1); 34 md.initialization.pressure = zeros(md.mesh.numberofvertices,1); 35 35 36 36 %Materials 37 md.initialization.temperature = (273.15-5.)*ones(md.mesh.numberofvertices,1) 38 md.initialization.waterfraction = zeros(md.mesh.numberofvertices,1) 39 md.initialization.watercolumn = zeros(md.mesh.numberofvertices,1) 40 md.materials.rheology_B = paterson(md.initialization.temperature) 41 md.materials.rheology_n = 3.*ones(md.mesh.numberofelements,1) 37 md.initialization.temperature = (273.15-5.)*ones(md.mesh.numberofvertices,1); 38 md.initialization.waterfraction = zeros(md.mesh.numberofvertices,1); 39 md.initialization.watercolumn = zeros(md.mesh.numberofvertices,1); 40 md.materials.rheology_B = paterson(md.initialization.temperature); 41 md.materials.rheology_n = 3.*ones(md.mesh.numberofelements,1); 42 42 43 43 %Thermal 44 md.thermal.isenthalpy = 0 45 md.thermal.spctemperature = NaN(md.mesh.numberofvertices,1) 44 md.thermal.isenthalpy = 0; 45 md.thermal.spctemperature = NaN(md.mesh.numberofvertices,1); 46 46 47 47 %Groundingline 48 md.groundingline.migration = 'SubelementMigration' 48 md.groundingline.migration = 'SubelementMigration'; 49 49 50 50 %Damage 51 md.damage.D = zeros(md.mesh.numberofvertices,1) 51 md.damage.D = zeros(md.mesh.numberofvertices,1); 52 52 53 53 %Surface mass balance and basal melting 54 md.surfaceforcings.mass_balance = 0.3*ones(md.mesh.numberofvertices,1) 55 md.basalforcings.melting_rate = md.surfaceforcings.mass_balance 54 md.surfaceforcings.mass_balance = 0.3*ones(md.mesh.numberofvertices,1); 55 md.basalforcings.melting_rate = md.surfaceforcings.mass_balance; 56 56 57 57 %Friction 58 md.friction.coefficient = 20.*ones(md.mesh.numberofvertices,1) 59 md.friction.coefficient(find(md.mask.groundedice_levelset<0.)) = 0. 60 md.friction.p = ones(md.mesh.numberofelements,1) 61 md.friction.q = ones(md.mesh.numberofelements,1) 58 md.friction.coefficient = 20.*ones(md.mesh.numberofvertices,1); 59 md.friction.coefficient(find(md.mask.groundedice_levelset<0.)) = 0.; 60 md.friction.p = ones(md.mesh.numberofelements,1); 61 md.friction.q = ones(md.mesh.numberofelements,1); 62 62 63 63 %Transient 64 md.transient.isstressbalance = 1 65 md.transient.islevelset = 1 66 md.transient.ismasstransport = 0 67 md.transient.isthermal = 0 68 md.transient.isgroundingline = 1 69 md.transient.isgia = 0 64 md.transient.isstressbalance = 1; 65 md.transient.islevelset = 1; 66 md.transient.ismasstransport = 0; 67 md.transient.isthermal = 0; 68 md.transient.isgroundingline = 1; 69 md.transient.isgia = 0; 70 70 71 71 %Stressbalance 72 md.stressbalance.maxiter = 100 73 md.stressbalance.viscosity_overshoot = 0.0 74 md.stressbalance.restol = 0.05 75 md.stressbalance.reltol = 0.05 76 md.stressbalance.abstol = NaN 72 md.stressbalance.maxiter = 100; 73 md.stressbalance.viscosity_overshoot = 0.0; 74 md.stressbalance.restol = 0.05; 75 md.stressbalance.reltol = 0.05; 76 md.stressbalance.abstol = NaN; 77 77 78 %Masstransport 79 md.masstransport.calvingrate = 0.*ones(md.mesh.numberofvertices,1) 80 md.masstransport.stabilization = 1. 78 %Masstransport; 79 md.masstransport.calvingrate = 0.*ones(md.mesh.numberofvertices,1); 80 md.masstransport.stabilization = 1.; 81 81 82 82 %Numerical parameters 83 md.thermal.stabilization = 1. 84 md.settings.waitonlock = 30 85 md.steadystate.reltol = 0.05 86 md.timestepping.time_step = 1. 87 md.timestepping.final_time = 3. 83 md.thermal.stabilization = 1.; 84 md.settings.waitonlock = 30; 85 md.steadystate.reltol = 0.05; 86 md.timestepping.time_step = 1.; 87 md.timestepping.final_time = 3.; 88 88 89 89 %Verbose 90 md.verbose = verbose(0) 90 md.verbose = verbose(0); 91 91 92 92 %Deal with boundary conditions: 93 md = SetIceShelfBC(md) 93 md = SetIceShelfBC(md); 94 94 95 %Change name so that no test have the same name 95 %Change name so that no test have the same name; 96 96 A = dbstack; 97 97 if (length(A)>2), md.miscellaneous.name=A(3).file(1:end-2); end
Note:
See TracChangeset
for help on using the changeset viewer.