Index: /issm/trunk-jpl/examples/AMR/mismip.py =================================================================== --- /issm/trunk-jpl/examples/AMR/mismip.py (revision 26558) +++ /issm/trunk-jpl/examples/AMR/mismip.py (revision 26558) @@ -0,0 +1,54 @@ +import numpy as np +from SetIceShelfBC import SetIceShelfBC + +# creating thickness +md.geometry.bed = -100 - np.abs(md.mesh.x) / 1000 +md.geometry.base = -90 * np.ones((md.mesh.numberofvertices)) +md.geometry.surface = 10 * np.ones((md.mesh.numberofvertices)) +md.geometry.thickness = md.geometry.surface - md.geometry.base +md.mask.ocean_levelset = -1 * np.ones((md.mesh.numberofvertices)) + +# creating basal drag +md.friction.coefficient = np.sqrt(10**7) * np.ones((md.mesh.numberofvertices)) #q = 1. +md.friction.p = 3 * np.ones((md.mesh.numberofelements)) +md.friction.q = np.zeros((md.mesh.numberofelements)) + +# creating flow law paramter +md.materials.rheology_B = 1 / ((1e-25)**(1 / 3)) * np.ones((md.mesh.numberofvertices)) +md.materials.rheology_n = 3 * np.ones((md.mesh.numberofelements)) +md.materials.rheology_law = 'None' + +# creating boundary conditions +md = SetIceShelfBC(md, './front.exp') +md.stressbalance.spcvx = np.nan * np.ones((md.mesh.numberofvertices)) +md.stressbalance.spcvy = np.nan * np.ones((md.mesh.numberofvertices)) +md.stressbalance.spcvz = np.nan * np.ones((md.mesh.numberofvertices)) +pos = np.nonzero(np.logical_or((np.abs(md.mesh.y - 50000) < 0.1), np.abs(md.mesh.y) < 0.1)) +md.stressbalance.spcvy[pos] = 0 +pos2 = np.nonzero(np.abs(md.mesh.x) < 0.1) +md.stressbalance.spcvx[pos2] = 0 +md.stressbalance.spcvz[pos] = np.nan +md.stressbalance.spcvz[pos2] = np.nan + +# creating forcing conditions +md.smb.mass_balance = 0.5 * np.ones((md.mesh.numberofvertices)) +md.basalforcings.geothermalflux = 0.5 * np.ones((md.mesh.numberofvertices)) +md.thermal.spctemperature = np.nan * np.ones((md.mesh.numberofvertices)) +md.groundingline.migration = 'SubelementMigration' + +# setting parameters +md.materials.rho_ice = 900 +md.materials.rho_water = 1000 +md.constants.g = 9.8 +md.constants.yts = 3600 * 24 * 365 +md.transient.isthermal = 0 +md.transient.isgroundingline = 1 +md.stressbalance.isnewton = 0 + +# setting inital condition +md.initialization.vx = np.ones((md.mesh.numberofvertices)) +md.initialization.vy = np.ones((md.mesh.numberofvertices)) +md.initialization.vz = np.ones((md.mesh.numberofvertices)) +md.initialization.vel = np.sqrt(2) * np.ones((md.mesh.numberofvertices)) +md.initialization.pressure = md.constants.g * md.materials.rho_ice * md.geometry.thickness +md.initialization.temperature = 273 * np.ones((md.mesh.numberofvertices)) Index: /issm/trunk-jpl/examples/AMR/runme.py =================================================================== --- /issm/trunk-jpl/examples/AMR/runme.py (revision 26558) +++ /issm/trunk-jpl/examples/AMR/runme.py (revision 26558) @@ -0,0 +1,78 @@ +# Mismip3D experiment with AMR using BAMG +from bamg import bamg +from model import * +from export_netCDF import export_netCDF +from setmask import setmask +from loadmodel import loadmodel +from parameterize import parameterize +from generic import generic +from exportVTK import exportVTK +from setflowequation import setflowequation +from solve import solve +from plotmodel import plotmodel + + +steps = [4] + +if 1 in steps: + print(' Step 1: Coarse mesh') + + #Generate an unstructured coarse mesh on the MISMIP domain with typical element edge length equal to 10, 000 m + md = bamg(model(), 'domain', './domain.exp', 'hmax', 10000, 'splitcorners', 1) + + export_netCDF(md, 'AMRCoarseMesh.nc') + +if 2 in steps: + print(' Step 2: Parameterization') + + md = loadmodel('AMRCoarseMesh.nc') + + md = setmask(md, '', '') + + # Run parameterization script to set up geometry, inital velocity, material properties, etc. + md = parameterize(md, './mismip.py') + + # Set the AMR properties and the refinement criteria + # Here, we are refining around the grounding line + # We impose the element resolution at grounding equal to 1000 m (1 km) + # The criterion used is the element distance to the grounding line + # The distance used here is 10000 m (10 km), used in both side around the grouding line (upstream and downstream) + md.amr.groundingline_resolution = 1000 + md.amr.groundingline_distance = 10000 + md.amr.hmin = 1000 # the same resolution used around the grounding line + md.amr.hmax = 10000 # the same coase resolution used to generate the coarse mesh + md.amr.gradation = 1.7 # this controls the ratio between two consecutive edges + md.amr.fieldname = 'None' # no field used here + md.amr.keepmetric = 0 # no field, no metric + + export_netCDF(md, 'AMRParam.nc') + +if 3 in steps: + print(' Step 3: Solve!') + + md = loadmodel('AMRParam.nc') + + # Run transient with adaptive mesh refinement + md.timestepping.time_step = 1 + md.timestepping.final_time = 500 # here, as example, only 500 yr. + md.settings.output_frequency = 10 # here, save results every 10 yr + md.stressbalance.maxiter = 30 + md.stressbalance.abstol = np.nan + md.stressbalance.restol = 1 + md.settings.solver_residue_threshold = 1e-2 # relaxing (the first stress balance solver iteration presents values higher than the original threshold. This probably happens because the initial velocity is set to one). + md.verbose = verbose('convergence', False, 'solution', True) + + # Specify that you want to run the model on your current (local host) computer + # Change the number of processors according to your machine (here np = 2) + md.cluster = generic('np', 2) + + # Set the AMR frequency, i.e., can be 1 or larger depending on how often the mesh needs to be updated + md.transient.amr_frequency = 1 # here, we are refining the mesh in every time step + + # Set the flow equation (SSA) and run + md = setflowequation(md, 'SSA', 'all') + md = solve(md, 'Transient') + + # Print the solutions and the mesh in VTK format (needs ParaView: https: / / www.paraview.org) + exportVTK('./VTKpy', md) + export_netCDF(md, 'AMRTransient.nc')