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