[issm-support] ISSM Questions
Morlighem, Mathieu (334H-Affiliate)
Mathieu.Morlighem at jpl.nasa.gov
Mon Feb 4 08:56:43 PST 2013
Hi Gordan,
1. The friction cannot be specified for our SIA model, the basal velocity is directly computed depending on the driving stress
Penta.cpp line 7528
constant_part=-1.58*pow((IssmDouble)10.0,-(IssmDouble)10.0)*rho_ice*gravity*thickness;
ub=constant_part*slope[0];
vb=constant_part*slope[1];
2. No, the Enthalpy is not coupled to a hydrology model. Whatever melts is directly added to the basal melting rates.
3. It looks like you have an old version of the code because vertex->dof does not exist anymore (it is now "pid" for parallel ID), but otherwise your code looks great !
If you want to update your code, go to trunk and "svn update", then you can reconfigure (http://issm.jpl.nasa.gov/download/unix/ last paragraph: ISSM Compilation) let us know if you encounter any problem.
You might want to use "TransientInputs" instead:
right now md.surfaceforcings.mass_balance is of size [md.mesh.numberofvertices 1], but you can use a time series as follows:
md.surfaceforcings.mass_balance =
[ smb(vertex1,time1) smb(vertex1,time2) …. smb(vertex1,timeN)
smb(vertex2,time2) smb(vertex2,time2) …. smb(vertex2,timeN)
…
smb(vertexn,time1) smb(vertexn,time2) …. smb(vertexn,timeN)
time1 tim2 …. timeN]
where time is in years. So in your case, if smb is your initial vector:
md.surfaceforcings.mass_balance = [smb smb+2; 0 2];
What will happen in the code is that:
- for time=0 it will use smb
- 0 < time < 2 yr it will use a linear interpolation between smb and smb+2 m/yr
- time >2 yr it will use smb+2 m/yr
Do not hesitate if you have any other question
Best,
Mathieu
--
Mathieu Morlighem
Assistant Project Scientist
Department of Earth System Science
University of California, Irvine
Croul Hall, Irvine, CA 92697-3100
(949) 824-1353 (UCI office)
(818) 354-4134 (JPL office)
(626) 429-5780 (cell)
Mathieu.Morlighem at uci.edu<mailto:Mathieu.Morlighem at uci.edu>
On Feb 3, 2013, at 10:54 AM, Gordan Stuhne <gordan at atmosp.physics.utoronto.ca<mailto:gordan at atmosp.physics.utoronto.ca>> wrote:
Hi Mathieu,
I have done some more investigation based on your
responses to my previous questions. The aphysical
melting rates do seem to be accompanied by numerical
noise, when I use the SSA-based inverse method
to specify basal friction in an SIA-based 3-D model. The
basal melting is still noisy in a HO-based 3-D model,
but the amplitude is more reasonable (O(10 cm/yr) max).
Presumably this is because the HO model includes
the SSA stress tensor components while the SIA model
does not.
I also tried your suggestion of using enthalpy
dynamics instead of temperature dynamics, but am
still working on the interpretation of the results.
Has the enthalpy model been fully coupled to
the hydrological model along the lines suggested
by Aschwanden et al (J. Glaciology, v. 58, 2012)?
On the technical side, is there some simple way
of exactly restarting simulations from disk
files? The only solution I can see is to create
a new initialization file by copying
md.model.TransientSolution... fields to
the appropriate places. Also, I would like
to use C++ code directly to change BCs at time
t inside the loop in the transient_core routine.
Would adding a call to the routine below
consistently work to, say, increase surface
mass balance by 2 m/a?
Best regards,
Gordan Stuhne.
--------------TEST CODE---------------------------
void my_addition(FemModel* femmodel){
double yts,*smass = NULL;
femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
int
nlocvert = femmodel->vertices->Size(),
nvert = femmodel->vertices->NumberOfVertices();
GetVectorFromInputsx(smass,femmodel->elements,
femmodel->nodes, femmodel->vertices,
femmodel->loads, femmodel->materials,
femmodel->parameters,
SurfaceforcingsMassBalanceEnum,
VertexEnum);
Vector *vsmass = new Vector(nvert);
for (int i=0;i<nlocvert;i++){
Vertex *vertex = (Vertex *)
femmodel->vertices->GetObjectByOffset(i);
int ind = vertex->dof;
vsmass->SetValue(ind,(*smass)[ind]+2.0/yts,INS_VAL);
}
vsmass->Assemble();
InputUpdateFromVectorx(femmodel->elements,
femmodel->nodes, femmodel->vertices,
femmodel->loads, femmodel->materials,
femmodel->parameters,vsmass,
SurfaceforcingsMassBalanceEnum,
VertexEnum);
xfree((void**)&smass);
xdelete(&vsmass);
}
------------------------------------------------------------------
On Mon, 2013-01-21 at 20:28 +0000, Morlighem, Mathieu (334H-Affiliate)
wrote:
Hello Gordan,
question 1:
We use the penalty method to constrain the ice temperature below
pressure melting point (Tpmp).
If the ice exceeds Tpmp, we add a penalty to the stiffness matrix and
load vector:
kappa x T = kappa x Tpmp
where kappa = kmax x 10^penalty_factor is a large number.
This forces T to be equal to Tpmp.
We can then retrieve the melting rate by computing how much energy was
necessary to constrain T:
kappa (T-Tpmp) = L / c Mb
where L is the latent heat, c is the heat capacity and Mb the melting
rate.
The problem with this method is that the melting rate is VERY noisy
and can reach values that are not physically realistic. This is due to
the fact that no regularization is applied (Mb is not in H1 but L2, if
you are familiar with Hilbert and Lebesgue spaces).
One alternative is to use the enthalpy model, which does not use
penalties. It is still under development but you could give it a try.
Just set md.thermal.isenthalpy=1;
Let us know if it works better.
question 2:
FS is much more sensitive than SIA numerically, and requires to have a
mesh that is refined enough, boundary conditions that are consistent
etc.
What we do generally is we don't jump from SIA to FS. You should try
SSA (MacAyeal) and HO (Blatter/Pattyn) first. Also, most of the time,
HO does a great job and it is not necessary to use a FS model.
Do not hesitate if you have more questions
Regards,
Mathieu
--
Mathieu Morlighem
Assistant Project Scientist
Department of Earth System Science
University of California, Irvine
Croul Hall, Irvine, CA 92697-3100
(949) 824-1353 (UCI office)
(818) 354-4134 (JPL office)
(626) 429-5780 (cell)
Mathieu.Morlighem at uci.edu<mailto:Mathieu.Morlighem at uci.edu>
On Jan 21, 2013, at 10:52 AM, Gordan Stuhne
<gordan at atmosp.physics.utoronto.ca> wrote:
Hello,
I am quite interested in your model, and have been working with it
using some of my own grids and data relating to Greenland.
My SIA (Hutter) results basically look rational for 100yr prognostic
runs (using thermal penalty_lock=10 for stability), but I have a few
questions. Firstly, the basal melting rates reported in the
output files are very high (on the order of m/yr). This is
presumably computed from latent heat fluxes where the
basal ice becomes temperate, but the details don't seem
to be documented in the code and I am not sure if I am
interpreting it right.
Secondly, I am having trouble trying to run the same case in
FS mode. I've tried adjusting time-step, viscosity overshoot,
and other parameters with no success. Is there some
useful rule of thumb for preparing the initial conditions?
Your insight would be greatly appreciated,
Thanks in advance,
Gordan Stuhne,
Dept. of Physics,
University of Toronto
_______________________________________________
issm-support mailing list
issm-support at issm.ess.uci.edu
http://issm.ess.uci.edu/cgi-bin/mailman/listinfo/issm-support
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://issm.ess.uci.edu/pipermail/issm-support/attachments/20130204/b9deb313/attachment.html>
More information about the issm-support
mailing list