Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 26254) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 26255) @@ -10,4 +10,6 @@ #include "../modules/modules.h" #include "../solutionsequences/solutionsequences.h" +#include "../classes/Inputs/DatasetInput.h" + void LevelsetAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ @@ -49,4 +51,5 @@ iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); +// iomodel->FetchDataToInput(inputs,elements,"md.mask.ref_levelset",MaskIceRefLevelsetEnum); iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum); iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyEnum); @@ -67,4 +70,10 @@ iomodel->FetchDataToInput(inputs,elements,"md.geometry.bed",BedEnum); break; + case CalvingTestEnum: + iomodel->FetchDataToInput(inputs,elements,"md.calving.stress_threshold_groundedice",CalvingStressThresholdGroundediceEnum); + iomodel->FetchDataToInput(inputs,elements,"md.calving.stress_threshold_floatingice",CalvingStressThresholdFloatingiceEnum); + iomodel->FetchDataToInput(inputs,elements,"md.calving.Qrates",CalvingNonlinearQRateEnum); + iomodel->FetchDataToInput(inputs,elements,"md.geometry.bed",BedEnum); + break; case CalvingMinthicknessEnum: iomodel->FetchDataToInput(inputs,elements,"md.geometry.bed",BedEnum); @@ -103,4 +112,7 @@ void LevelsetAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ + IssmDouble* coefficients=NULL; + int nco; + parameters->AddObject(iomodel->CopyConstantObject("md.levelset.stabilization",LevelsetStabilizationEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.levelset.reinit_frequency",LevelsetReinitFrequencyEnum)); @@ -117,4 +129,13 @@ parameters->AddObject(iomodel->CopyConstantObject("md.calving.min_thickness",CalvingMinthicknessEnum)); break; + case CalvingTestEnum: + parameters->AddObject(iomodel->CopyConstantObject("md.calving.min_thickness",CalvingMinthicknessEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.calving.dt_order",CalvingDtimeorderEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.calving.nonlinear_law",CalvingNonlinearLawEnum)); + + /* Add the coefficients for the new nonlinear law */ + iomodel->FetchData(&coefficients, &nco, NULL, "md.calving.coefficients"); + parameters->AddObject(new DoubleMatParam(CalvingNonlinearCoeffEnum, coefficients, nco, 1)); + xDelete(coefficients); case CalvingMinthicknessEnum: parameters->AddObject(iomodel->CopyConstantObject("md.calving.min_thickness",CalvingMinthicknessEnum)); @@ -215,4 +236,5 @@ IssmDouble* dbasis = xNew(2*numnodes); IssmDouble* Bprime = NULL; + IssmDouble tau; if(stabilization==2){ Bprime = xNew(dim*numnodes); @@ -340,8 +362,9 @@ /*Intermediaries */ - int domaintype; + int domaintype, stabilization; IssmDouble Jdet,dt; IssmDouble lsf; IssmDouble* xyz_list = NULL; + IssmDouble vx,vy,vel,h,tau; /*Fetch number of nodes and dof for this finite element*/ @@ -351,4 +374,5 @@ ElementVector* pe = basalelement->NewElementVector(); basalelement->FindParam(&dt,TimesteppingTimeStepEnum); + basalelement->FindParam(&stabilization,LevelsetStabilizationEnum); if(dt!=0.){ @@ -359,4 +383,6 @@ basalelement->GetVerticesCoordinates(&xyz_list); Input* levelset_input = basalelement->GetInput(MaskIceLevelsetEnum); _assert_(levelset_input); + + // h=element->CharacteristicLength(); /* Start looping on the number of gaussian points: */ @@ -473,4 +499,39 @@ } + if(calvinglaw==CalvingTestEnum){ + + /*Get minimum thickness threshold*/ + femmodel->parameters->FindParam(&min_thickness,CalvingMinthicknessEnum); + + /*Loop over all elements of this partition*/ + for(Object* & object : femmodel->elements->objects){ + Element* element = xDynamicCast(object); + + int numnodes = element->GetNumberOfNodes(); + Gauss* gauss = element->NewGauss(); + Input* H_input = element->GetInput(ThicknessEnum); _assert_(H_input); + Input* b_input = element->GetInput(BedEnum); _assert_(b_input); + Input* sl_input = element->GetInput(SealevelEnum); _assert_(sl_input); + + /*Potentially constrain nodes of this element*/ + for(int in=0;inGaussNode(element->GetElementType(),in); + Node* node=element->GetNode(in); + if(!node->IsActive()) continue; + + H_input->GetInputValue(&thickness,gauss); + b_input->GetInputValue(&bed,gauss); + sl_input->GetInputValue(&sealevel,gauss); + if(thicknessApplyConstraint(0,+1.); + } + else { + /* no ice, set no spc */ + node->DofInFSet(0); + } + } + delete gauss; + } + } if(calvinglaw==CalvingHabEnum){ Index: /issm/trunk-jpl/src/c/classes/Cflevelsetmisfit.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Cflevelsetmisfit.cpp (revision 26254) +++ /issm/trunk-jpl/src/c/classes/Cflevelsetmisfit.cpp (revision 26255) @@ -89,9 +89,7 @@ marshallhandle->call(this->definitionenum); - marshallhandle->call(this->local); marshallhandle->call(this->model_enum); marshallhandle->call(this->name); marshallhandle->call(this->observation_enum); - marshallhandle->call(this->timeinterpolation); marshallhandle->call(this->weights_enum); marshallhandle->call(this->datatime); Index: /issm/trunk-jpl/src/c/datastructures/DataSet.cpp =================================================================== --- /issm/trunk-jpl/src/c/datastructures/DataSet.cpp (revision 26254) +++ /issm/trunk-jpl/src/c/datastructures/DataSet.cpp (revision 26255) @@ -251,4 +251,9 @@ this->AddObject(res); } + else if(obj_enum==CflevelsetmisfitEnum){ + Cflevelsetmisfit* Cflevelset=new Cflevelsetmisfit(); + Cflevelset->Marshall(marshallhandle); + this->AddObject(Cflevelset); + } else _error_("could not recognize enum type: " << obj_enum << ": " << EnumToStringx(obj_enum) ); } Index: /issm/trunk-jpl/test/NightlyRun/test3203.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test3203.m (revision 26255) +++ /issm/trunk-jpl/test/NightlyRun/test3203.m (revision 26255) @@ -0,0 +1,91 @@ +%Test Name: SquareShelfTransientLevelsetMisfitcodipack + +md=triangle(model(),'../Exp/Square.exp',50000.); +md=setmask(md,'all',''); +md=parameterize(md,'../Par/SquareShelf.par'); +md=setflowequation(md,'SSA','all'); +md.cluster=generic('name',oshostname(),'np',3); + +%Do not kill ice bergs as all is floating +md.levelset.kill_icebergs=0; + +x = md.mesh.x; +xmin = min(x); +xmax = max(x); +Lx = (xmax-xmin); +alpha = 2./3.; +md.mask.ice_levelset = ((x - alpha*Lx)>0) - ((x - alpha*Lx)<0); + +md.timestepping.time_step=10; +md.timestepping.final_time=30; + +%Transient +md.transient.isstressbalance=1; +md.transient.ismasstransport=1; +md.transient.issmb=1; +md.transient.isthermal=0; +md.transient.isgroundingline=0; +md.transient.ismovingfront=1; + +md.calving=calvinglevermann(); +md.calving.coeff=4.89e13*ones(md.mesh.numberofvertices,1); +md.frontalforcings.meltingrate=zeros(md.mesh.numberofvertices,1); +md.levelset.spclevelset=NaN(md.mesh.numberofvertices,1); +md.levelset.migration_max = 1e8; + +md = solve(md,'tr'); +%plotmodel(md,'axis#all','tight','data',md.materials.rheology_B(1:end-1,1),'caxis#all',[ 1.3 1.9]*10^8,'title','"True" B',... +%'data',md.materials.rheology_B(1:end-1,2),'title','"True" B 2') + +%Modify rheology, now constant +md.materials.rheology_B(1:end-1,:) = 1.8e8; + +%Set cost function +weights= ones(md.mesh.numberofvertices,1); +count = 1; + +for i=1:numel(md.results.TransientSolution) + time = md.results.TransientSolution(i).time; + md.outputdefinition.definitions{count}=cflevelsetmisfit('name',['LevelsetMisfit' num2str(count)],... + 'definitionstring',['Outputdefinition' num2str(count)],... + 'model_string','MaskIceLevelset','observation_string','LevelsetObservation',... + 'observation',reinitializelevelset(md, md.results.TransientSolution(i).MaskIceLevelset),'weights',weights,'weights_string','WeightsLevelsetObservation',... + 'datatime',time); + md.autodiff.dependents{count} = dependent('name',['Outputdefinition' num2str(count)],'type','scalar','fos_reverse_index',1); + + count = count+1; +end + +%Independent +min_params = md.materials.rheology_B; min_params(1:end-1,:) = cuffey(273); +max_params = md.materials.rheology_B; max_params(1:end-1,:) = cuffey(200); +md.autodiff.independents{1} = independent('name','MaterialsRheologyBbar',... + 'md_name','md.materials.rheology_B',... + 'control_size',size(md.materials.rheology_B,2),... + 'type','vertex',... %Really needed?? + 'min_parameters',min_params,... + 'max_parameters',max_params,... + 'control_scaling_factor',1e8); + +md.inversion=adm1qn3inversion(md.inversion); +md.inversion.iscontrol=1; +md.inversion.maxiter=4; +md.inversion.maxsteps=md.inversion.maxiter; +md.inversion.dxmin=1e-5; +md.autodiff.isautodiff=1; +md.autodiff.driver='fos_reverse'; + +%Go solve! +md.verbose=verbose(0); +md=solve(md,'tr'); +%plotmodel(md,'axis#all','tight','data',md.results.TransientSolution(1).MaterialsRheologyBbar(:,1),'caxis#all',[ 1.3 1.9]*10^8,'title','B1',... +%'data',md.results.TransientSolution(1).MaterialsRheologyBbar(:,2),'title','B2') + +%Fields and tolerances to track changes +field_names ={'Gradient','Misfit','Rheology'}; +field_tolerances={1e-12,1e-12,1e-12}; +field_values={... + (md.results.TransientSolution(1).Gradient1),... + (md.results.TransientSolution(1).J),... + (md.results.TransientSolution(1).MaterialsRheologyBbar),... + };