Index: /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp (revision 26893) +++ /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp (revision 26894) @@ -73,5 +73,5 @@ } - iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum); + iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum); iomodel->FetchDataToInput(inputs,elements,"md.initialization.sealevel",SealevelEnum,0); iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); @@ -96,7 +96,7 @@ iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.floatingice_melting_rate",BasalforcingsFloatingiceMeltingRateEnum); if(isstochastic){ - iomodel->FetchDataToInput(inputs,elements,"md.stochasticforcing.default_id",StochasticForcingDefaultIdEnum); - iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.floatingice_melting_rate",BaselineBasalforcingsFloatingiceMeltingRateEnum); - } + iomodel->FetchDataToInput(inputs,elements,"md.stochasticforcing.default_id",StochasticForcingDefaultIdEnum); + iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.floatingice_melting_rate",BaselineBasalforcingsFloatingiceMeltingRateEnum); + } break; case LinearFloatingMeltRateEnum: @@ -113,7 +113,7 @@ iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.upperwater_elevation",BasalforcingsSpatialUpperwaterElevationEnum); if(isstochastic){ - iomodel->FetchDataToInput(inputs,elements,"md.stochasticforcing.default_id",StochasticForcingDefaultIdEnum); - iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_melting_rate",BaselineBasalforcingsSpatialDeepwaterMeltingRateEnum); - } + iomodel->FetchDataToInput(inputs,elements,"md.stochasticforcing.default_id",StochasticForcingDefaultIdEnum); + iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_melting_rate",BaselineBasalforcingsSpatialDeepwaterMeltingRateEnum); + } break; case BasalforcingsPicoEnum: @@ -148,5 +148,5 @@ }/*}}}*/ ElementMatrix* FreeSurfaceBaseAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/ -_error_("Not implemented"); + _error_("Not implemented"); }/*}}}*/ ElementMatrix* FreeSurfaceBaseAnalysis::CreateKMatrix(Element* element){/*{{{*/ @@ -157,5 +157,5 @@ IssmDouble *xyz_list = NULL; IssmDouble Jdet,D_scalar,dt,h; - IssmDouble vel,vx,vy; + IssmDouble vel,vx,vy,tau; /*Get basal element*/ @@ -230,19 +230,5 @@ } - if(stabilization==2){ - /*Streamline upwinding*/ - if(dim==1){ - vel=fabs(vx)+1.e-8; - D[0] = h/(2.*vel)*vx*vx; - } - else{ - vel=sqrt(vx*vx+vy*vy)+1.e-8; - D[0*dim+0]=h/(2*vel)*vx*vx; - D[1*dim+0]=h/(2*vel)*vy*vx; - D[0*dim+1]=h/(2*vel)*vx*vy; - D[1*dim+1]=h/(2*vel)*vy*vy; - } - } - else if(stabilization==1){ + if(stabilization==1){ /*SSA*/ if(dim==1){ @@ -255,4 +241,31 @@ D[0*dim+0]=h/2.0*fabs(vx); D[1*dim+1]=h/2.0*fabs(vy); + } + } + else if(stabilization==2){ + /*Streamline upwinding*/ + if(dim==1){ + vel=fabs(vx)+1.e-8; + D[0] = h/(2.*vel)*vx*vx; + } + else{ + vel=sqrt(vx*vx+vy*vy)+1.e-8; + D[0*dim+0]=h/(2*vel)*vx*vx; + D[1*dim+0]=h/(2*vel)*vy*vx; + D[0*dim+1]=h/(2*vel)*vx*vy; + D[1*dim+1]=h/(2*vel)*vy*vy; + } + + } + else if(stabilization==5){ + /*SUPG*/ + if(dim==1){ + vx_input->GetInputAverage(&vx); + tau=h/(2.*fabs(vx)+1e-10); + } + else{ + vx_input->GetInputAverage(&vx); + vy_input->GetInputAverage(&vy); + tau=1*h/(2.*pow(vx*vx+vy*vy,0.5)+1e-10); } } @@ -263,7 +276,7 @@ for(int j=0;jvalues[i*numnodes+j] += ( - dbasis[0*numnodes+i] *(D[0*dim+0]*dbasis[0*numnodes+j] + D[0*dim+1]*dbasis[1*numnodes+j]) + - dbasis[1*numnodes+i] *(D[1*dim+0]*dbasis[0*numnodes+j] + D[1*dim+1]*dbasis[1*numnodes+j]) - ); + dbasis[0*numnodes+i] *(D[0*dim+0]*dbasis[0*numnodes+j] + D[0*dim+1]*dbasis[1*numnodes+j]) + + dbasis[1*numnodes+i] *(D[1*dim+0]*dbasis[0*numnodes+j] + D[1*dim+1]*dbasis[1*numnodes+j]) + ); } } @@ -271,4 +284,19 @@ else{ for(int i=0;ivalues[i*numnodes+j] += dbasis[0*numnodes+i]*D[0]*dbasis[0*numnodes+j]; + } + } + else if(stabilization==5){ + D_scalar=gauss->weight*Jdet*dt; + if(dim==2){ + for(int i=0;ivalues[i*numnodes+j]+=tau*D_scalar* + (vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i])* + (vx*dbasis[0*numnodes+j]+vy*dbasis[1*numnodes+j]); + } + } + } + else{ + for(int i=0;ivalues[i*numnodes+j]+=tau*D_scalar*(vx*dbasis[0*numnodes+i])*(vx*dbasis[0*numnodes+j]); } } @@ -285,8 +313,9 @@ }/*}}}*/ ElementVector* FreeSurfaceBaseAnalysis::CreatePVector(Element* element){/*{{{*/ + /*Intermediaries*/ - int domaintype,dim; + int domaintype,dim,stabilization; IssmDouble Jdet,dt; - IssmDouble gmb,fmb,mb,bed,phi,vz; + IssmDouble gmb,fmb,mb,bed,vx,vy,vz,tau; Element* basalelement = NULL; IssmDouble *xyz_list = NULL; @@ -314,8 +343,13 @@ /*Fetch number of nodes and dof for this finite element*/ int numnodes = basalelement->GetNumberOfNodes(); + int melt_style,point1; + IssmDouble fraction1,fraction2; + bool mainlyfloating; /*Initialize Element vector and other vectors*/ ElementVector* pe = basalelement->NewElementVector(); IssmDouble* basis = xNew(numnodes); + IssmDouble* dbasis = xNew(dim*numnodes); + IssmDouble gllevelset,phi=1.; /*Retrieve all inputs and parameters*/ @@ -326,13 +360,35 @@ Input* fmb_input = basalelement->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(fmb_input); Input* base_input = basalelement->GetInput(BaseEnum); _assert_(base_input); - Input* vz_input = NULL; + Input* gllevelset_input = basalelement->GetInput(MaskOceanLevelsetEnum); _assert_(gllevelset_input); + Input* vz_input = NULL; + Input* vx_input = NULL; + Input* vy_input = NULL; switch(dim){ - case 1: vz_input = basalelement->GetInput(VyEnum); _assert_(vz_input); break; - case 2: vz_input = basalelement->GetInput(VzEnum); _assert_(vz_input); break; + case 1: + vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input); + vz_input=basalelement->GetInput(VyEnum) ; _assert_(vz_input); + break; + case 2: + vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input); + vy_input=basalelement->GetInput(VyEnum); _assert_(vy_input); + vz_input=basalelement->GetInput(VzEnum); _assert_(vz_input); + break; default: _error_("not implemented"); } + IssmDouble h = basalelement->CharacteristicLength(); + + /*Recover portion of element that is grounded*/ + basalelement->GetVerticesCoordinates(&xyz_list); + phi=basalelement->GetGroundedPortion(xyz_list); + Gauss* gauss = NULL; + if(melt_style==SubelementMelt2Enum){ + basalelement->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating); + gauss = basalelement->NewGauss(point1,fraction1,fraction2,3); + } + else{ + gauss = basalelement->NewGauss(3); + } /* Start looping on the number of gaussian points: */ - Gauss* gauss=basalelement->NewGauss(2); while(gauss->next()){ @@ -340,13 +396,49 @@ basalelement->NodalFunctions(basis,gauss); - vz_input->GetInputValue(&vz,gauss); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + vz_input->GetInputValue(&vz,gauss); gmb_input->GetInputValue(&gmb,gauss); fmb_input->GetInputValue(&fmb,gauss); base_input->GetInputValue(&bed,gauss); groundedice_input->GetInputValue(&phi,gauss); - if(phi>0) mb=gmb; - else mb=fmb; + gllevelset_input->GetInputValue(&gllevelset,gauss); + if(melt_style==SubelementMelt1Enum){ + //if (phi>0.999999999) mb=gmb; + //else mb=(1-phi)*fmb+phi*gmb; // phi is the fraction of grounded ice so (1-phi) is floating + if(phi>0) mb=gmb; + else mb=fmb; + } + else if(melt_style==SubelementMelt2Enum){ + if(gllevelset>0.) mb=gmb; + else mb=fmb; + } + else if(melt_style==NoMeltOnPartiallyFloatingEnum){ + if (phi<0.00000001) mb=fmb; + else mb=gmb; + } + else if(melt_style==FullMeltOnPartiallyFloatingEnum){ + if (phi<0.99999999) mb=fmb; + else mb=gmb; + } + else _error_("melt interpolation "<values[i]+=Jdet*gauss->weight*(bed+dt*(mb) + dt*vz)*basis[i]; + + if(stabilization==5){ + /*SUPG*/ + basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + if(dim==1){ + vx_input->GetInputAverage(&vx); + tau=h/(2.*fabs(vx)+1e-10); + for(int i=0;ivalues[i]+=Jdet*gauss->weight*(dt*mb+dt*vz)*tau*(vx*dbasis[0*numnodes+i]); + } + else{ + vx_input->GetInputAverage(&vx); + vy_input->GetInputAverage(&vy); + tau=1*h/(2.*pow(vx*vx+vy*vy,0.5)+1e-10); + for(int i=0;ivalues[i]+=Jdet*gauss->weight*(bed*0.+dt*mb+dt*vz)*tau*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i]); + } + } } @@ -360,5 +452,5 @@ }/*}}}*/ void FreeSurfaceBaseAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ - _error_("not implemented yet"); + _error_("not implemented yet"); }/*}}}*/ void FreeSurfaceBaseAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_interp,int control_index){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp (revision 26893) +++ /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp (revision 26894) @@ -80,4 +80,5 @@ iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum); + iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyEnum); if(iomodel->domaintype!=Domain2DhorizontalEnum){ iomodel->FetchDataToInput(inputs,elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum); @@ -111,5 +112,5 @@ }/*}}}*/ ElementMatrix* FreeSurfaceTopAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/ -_error_("Not implemented"); + _error_("Not implemented"); }/*}}}*/ ElementMatrix* FreeSurfaceTopAnalysis::CreateKMatrix(Element* element){/*{{{*/ @@ -120,5 +121,5 @@ IssmDouble *xyz_list = NULL; IssmDouble Jdet,D_scalar,dt,h; - IssmDouble vel,vx,vy; + IssmDouble vel,vx,vy,tau; /*Get top element*/ @@ -195,5 +196,19 @@ } - if(stabilization==2){ + if(stabilization==1){ + /*artifical diffusion*/ + if(dim==1){ + vx_input->GetInputAverage(&vx); + D[0]=h/2.*fabs(vx); + } + else{ + vx_input->GetInputAverage(&vx); + vy_input->GetInputAverage(&vy); + + D[0*dim+0]=h/2.0*fabs(vx); + D[1*dim+1]=h/2.0*fabs(vy); + } + } + else if(stabilization==2){ /*Streamline upwinding*/ if(dim==1){ @@ -209,16 +224,14 @@ } } - else if(stabilization==1){ - /*SSA*/ + else if(stabilization==5){ + /*SUPG*/ if(dim==1){ vx_input->GetInputAverage(&vx); - D[0]=h/2.*fabs(vx); + tau=h/(2.*fabs(vx)+1e-10); } else{ vx_input->GetInputAverage(&vx); vy_input->GetInputAverage(&vy); - - D[0*dim+0]=h/2.0*fabs(vx); - D[1*dim+1]=h/2.0*fabs(vy); + tau=1*h/(2.*pow(vx*vx+vy*vy,0.5)+1e-10); } } @@ -229,7 +242,7 @@ for(int j=0;jvalues[i*numnodes+j] += ( - dbasis[0*numnodes+i] *(D[0*dim+0]*dbasis[0*numnodes+j] + D[0*dim+1]*dbasis[1*numnodes+j]) + - dbasis[1*numnodes+i] *(D[1*dim+0]*dbasis[0*numnodes+j] + D[1*dim+1]*dbasis[1*numnodes+j]) - ); + dbasis[0*numnodes+i] *(D[0*dim+0]*dbasis[0*numnodes+j] + D[0*dim+1]*dbasis[1*numnodes+j]) + + dbasis[1*numnodes+i] *(D[1*dim+0]*dbasis[0*numnodes+j] + D[1*dim+1]*dbasis[1*numnodes+j]) + ); } } @@ -237,4 +250,19 @@ else{ for(int i=0;ivalues[i*numnodes+j] += dbasis[0*numnodes+i]*D[0]*dbasis[0*numnodes+j]; + } + } + else if(stabilization==5){ + D_scalar=gauss->weight*Jdet*dt; + if(dim==2){ + for(int i=0;ivalues[i*numnodes+j]+=tau*D_scalar* + (vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i])* + (vx*dbasis[0*numnodes+j]+vy*dbasis[1*numnodes+j]); + } + } + } + else{ + for(int i=0;ivalues[i*numnodes+j]+=tau*D_scalar*(vx*dbasis[0*numnodes+i])*(vx*dbasis[0*numnodes+j]); } } @@ -251,8 +279,9 @@ }/*}}}*/ ElementVector* FreeSurfaceTopAnalysis::CreatePVector(Element* element){/*{{{*/ + /*Intermediaries*/ - int domaintype,dim; + int domaintype,dim,stabilization; IssmDouble Jdet,dt; - IssmDouble ms,surface,vz; + IssmDouble ms,surface,vx,vy,vz,tau; Element* topelement = NULL; IssmDouble *xyz_list = NULL; @@ -284,4 +313,5 @@ ElementVector* pe = topelement->NewElementVector(); IssmDouble* basis = xNew(numnodes); + IssmDouble* dbasis = xNew(dim*numnodes); /*Retrieve all inputs and parameters*/ @@ -291,9 +321,19 @@ Input* surface_input = topelement->GetInput(SurfaceEnum); _assert_(surface_input); Input* vz_input = NULL; + Input* vx_input = NULL; + Input* vy_input = NULL; switch(dim){ - case 1: vz_input = topelement->GetInput(VyEnum); _assert_(vz_input); break; - case 2: vz_input = topelement->GetInput(VzEnum); _assert_(vz_input); break; + case 1: + vx_input=topelement->GetInput(VxEnum); _assert_(vx_input); + vz_input = topelement->GetInput(VyEnum) ; _assert_(vz_input); + break; + case 2: + vx_input=topelement->GetInput(VxEnum); _assert_(vx_input); + vy_input = topelement->GetInput(VyEnum); _assert_(vy_input); + vz_input = topelement->GetInput(VzEnum); _assert_(vz_input); + break; default: _error_("not implemented"); } + IssmDouble h = topelement->CharacteristicLength(); /*Initialize mb_correction to 0, do not forget!:*/ @@ -312,7 +352,24 @@ } + if(stabilization==5){ + /*SUPG*/ + topelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + if(dim==1){ + vx_input->GetInputAverage(&vx); + tau=h/(2.*fabs(vx)+1e-10); + for(int i=0;ivalues[i]+=Jdet*gauss->weight*(dt*ms+dt*vz)*tau*(vx*dbasis[0*numnodes+i]); + } + else{ + vx_input->GetInputAverage(&vx); + vy_input->GetInputAverage(&vy); + tau=h/(2.*pow(vx*vx+vy*vy,0.5)+1e-10); + for(int i=0;ivalues[i]+=Jdet*gauss->weight*(dt*ms+dt*vz)*tau*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i]); + } + } + /*Clean up and return*/ xDelete(xyz_list); xDelete(basis); + xDelete(dbasis); delete gauss; if(topelement->IsSpawnedElement()){topelement->DeleteMaterials(); delete topelement;}; @@ -321,5 +378,5 @@ }/*}}}*/ void FreeSurfaceTopAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ - _error_("not implemented yet"); + _error_("not implemented yet"); }/*}}}*/ void FreeSurfaceTopAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_interp,int control_index){/*{{{*/ @@ -329,4 +386,71 @@ element->InputUpdateFromSolutionOneDof(solution,SurfaceEnum); + + /*Now, we need to do some "processing"*/ + int numvertices = element->GetNumberOfVertices(); + int migration_style; + + IssmDouble* surface = xNew(numvertices); + IssmDouble* newsurface = xNew(numvertices); + IssmDouble* thickness = xNew(numvertices); + IssmDouble* base = xNew(numvertices); + IssmDouble* bed = xNew(numvertices); + IssmDouble* phi = xNew(numvertices); + IssmDouble* sealevel = xNew(numvertices); + + IssmDouble minthickness = element->FindParam(MasstransportMinThicknessEnum); + IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum); + IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum); + + bool isgroundingline; + element->FindParam(&isgroundingline,TransientIsgroundinglineEnum); + element->FindParam(&migration_style,GroundinglineMigrationEnum); + if(isgroundingline) element->GetInputListOnVertices(&bed[0],BedEnum); + + element->GetInputListOnVertices(&base[0],BaseEnum); + element->GetInputListOnVertices(&surface[0],SurfaceEnum); + element->GetInputListOnVertices(&phi[0],MaskOceanLevelsetEnum); + element->GetInputListOnVertices(&sealevel[0],SealevelEnum); + + for(int i=0;i(thickness[i])) _error_("NaN found in solution vector"); + if(xIsInf(thickness[i])) _error_("Inf found in solution vector"); + + /* check for thickness0.){ + if(base[i]<=bed[i]) base[i] = bed[i]; + newsurface[i] = base[i]+minthickness; + }else{ + // assume floatation condition + newsurface[i] = (1.-rho_ice/rho_water)*minthickness; + base[i] = -rho_ice/rho_water*minthickness; + } + } + + /* update thickness */ + thickness[i]=newsurface[i]-base[i]; + /* some checks */ + if(thickness[i]<0.) _error_("thickness<0"); + if(newsurface[i]AddInput(BaseEnum,base,element->GetElementType()); + element->AddInput(SurfaceEnum,newsurface,element->GetElementType()); + element->AddInput(ThicknessEnum,thickness,element->GetElementType()); + + /* Free resources */ + xDelete(newsurface); + xDelete(surface); + xDelete(thickness); + xDelete(base); + xDelete(bed); + xDelete(phi); + xDelete(sealevel); }/*}}}*/ void FreeSurfaceTopAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ Index: /issm/trunk-jpl/src/c/cores/masstransport_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/masstransport_core.cpp (revision 26893) +++ /issm/trunk-jpl/src/c/cores/masstransport_core.cpp (revision 26894) @@ -49,4 +49,8 @@ solutionsequence_linear(femmodel); femmodel->parameters->SetParam(SurfaceEnum,InputToExtrudeEnum); + extrudefromtop_core(femmodel); + femmodel->parameters->SetParam(ThicknessEnum,InputToExtrudeEnum); + extrudefromtop_core(femmodel); + femmodel->parameters->SetParam(BaseEnum,InputToExtrudeEnum); extrudefromtop_core(femmodel); }