Index: /issm/trunk-jpl/src/c/cores/movingfront_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/movingfront_core.cpp (revision 26972) +++ /issm/trunk-jpl/src/c/cores/movingfront_core.cpp (revision 26973) @@ -16,5 +16,5 @@ /* intermediaries */ - bool save_results,isstressbalance,ismasstransport,isthermal,isenthalpy,islevelset,ismovingfront,killicebergs,haskillberg; + bool save_results,isstressbalance,ismasstransport,isthermal,isenthalpy,islevelset,ismovingfront,killicebergs; int domaintype, num_extrapol_vars, index,reinit_frequency,step; int* extrapol_vars=NULL; @@ -120,9 +120,14 @@ /*Kill ice berg to avoid free body motion*/ if(killicebergs){ + int killberg = 0; if(VerboseSolution()) _printf0_(" looking for icebergs to kill\n"); - KillIcebergsx(femmodel); - femmodel->parameters->FindParam(&haskillberg,LevelsetHasKillIcebergsEnum); - if (haskillberg) { - if(VerboseSolution()) _printf0_(" reinitializing level set after killing icebergs\n"); + killberg = KillIcebergsx(femmodel); + /*wait for all cores*/ + int totalkill; + ISSM_MPI_Reduce(&killberg,&totalkill,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&totalkill,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + + if (totalkill > 0) { + if(VerboseSolution()) _printf0_(" reinitializing level set after killing " << totalkill << " icebergs\n"); femmodel->ResetLevelset(); ResetBoundaryConditions(femmodel,LevelsetAnalysisEnum); Index: /issm/trunk-jpl/src/c/modules/KillIcebergsx/KillIcebergsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/KillIcebergsx/KillIcebergsx.cpp (revision 26972) +++ /issm/trunk-jpl/src/c/modules/KillIcebergsx/KillIcebergsx.cpp (revision 26973) @@ -9,5 +9,5 @@ #include "../InputUpdateFromVectorx/InputUpdateFromVectorx.h" -void KillIcebergsx(FemModel* femmodel){ +int KillIcebergsx(FemModel* femmodel){ /*Intermediaries*/ @@ -18,4 +18,5 @@ IssmDouble *local_mask = xNewZeroInit(nbv_local); bool *element_flag = xNewZeroInit(femmodel->elements->Size()); + IssmDouble ice; /*Step 1, go through all elements and put 1 in local_mask if the element is grounded*/ @@ -27,8 +28,18 @@ element_flag[i] = true; } - else{ + else { if(element->IsAllGrounded()){ + /* only look at element with ice but not fully grounded */ int numvertices = element->GetNumberOfVertices(); - for(int v=0;vvertices[v]->Lid()] = 1.; + Gauss* gauss=element->NewGauss(); + Input* icelevel_input = element->GetInput(MaskIceLevelsetEnum); _assert_(icelevel_input); + + for(int v=0;vGaussVertex(v); + icelevel_input->GetInputValue(&ice,gauss); + /* The initial mask is very strict, we look at all grounded elements and set the mask for ice nodes only. */ + if (ice < 0) local_mask[element->vertices[v]->Lid()] = 1.; + } + delete gauss; } } @@ -48,5 +59,4 @@ /*Local iterations on partition*/ bool keepgoing = true; - int didsomething = 0; int iter = 1; while(keepgoing){ @@ -54,5 +64,4 @@ keepgoing = false; - didsomething = 0; int i=0; for(Object* & object : femmodel->elements->objects){ @@ -62,7 +71,11 @@ int numvertices = element->GetNumberOfVertices(); bool found1 = false; + IssmDouble sumlocalmask = 0.; + for(int j=0;jvertices[j]->Lid(); - if(local_mask[lid]>0.){ + /*we need to have at least a sharing edge, to extend the mask*/ + sumlocalmask += local_mask[lid]; + if(sumlocalmask > 1.5){ found1 = true; break; @@ -73,10 +86,7 @@ for(int j=0;jvertices[j]->Lid(); - if(local_mask[lid]==0.){ - local_mask[lid]=1.; - keepgoing = true; - didsomething = 1; - } + local_mask[lid]=1.; } + keepgoing = true; } } @@ -99,5 +109,5 @@ xDelete(element_flag); - bool killbergReinit = false; + int killbergReinit = 0; /*OK, now deactivate iceberg and count the number of deactivated vertices*/ for(Object* & object : femmodel->elements->objects){ @@ -121,12 +131,12 @@ element->AddInput(MaskIceLevelsetEnum,values,P1Enum); xDelete(values); - killbergReinit = true; + killbergReinit += 1; } } } - /*Recompute the sign distance for the levelset function*/ - femmodel->parameters->SetParam(killbergReinit,LevelsetHasKillIcebergsEnum); - /*cleanup*/ xDelete(local_mask); + + /*Recompute the sign distance for the levelset function*/ + return killbergReinit; } Index: /issm/trunk-jpl/src/c/modules/KillIcebergsx/KillIcebergsx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/KillIcebergsx/KillIcebergsx.h (revision 26972) +++ /issm/trunk-jpl/src/c/modules/KillIcebergsx/KillIcebergsx.h (revision 26973) @@ -6,5 +6,5 @@ /* local prototypes: */ -void KillIcebergsx(FemModel* femmodel); +int KillIcebergsx(FemModel* femmodel); #endif Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26972) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26973) @@ -270,5 +270,4 @@ LevelsetReinitFrequencyEnum, LevelsetStabilizationEnum, - LevelsetHasKillIcebergsEnum, LockFileNameEnum, LoveAllowLayerDeletionEnum, Index: /issm/trunk-jpl/src/m/parameterization/killberg.m =================================================================== --- /issm/trunk-jpl/src/m/parameterization/killberg.m (revision 26972) +++ /issm/trunk-jpl/src/m/parameterization/killberg.m (revision 26973) @@ -28,6 +28,7 @@ %isgrounded = max(md.mask.ocean_levelset(md.mesh.elements),[],2)>0; pos = find(isgrounded); -element_flag(pos) = 1; +%element_flag(pos) = 1; mask(md.mesh.elements(pos,:)) = 1; +mask(md.mask.ice_levelset>=0) = 0; iter = 1; Index: /issm/trunk-jpl/test/NightlyRun/test813.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test813.m (revision 26973) +++ /issm/trunk-jpl/test/NightlyRun/test813.m (revision 26973) @@ -0,0 +1,59 @@ +%Test Name: SquareShelfLevelsetKillberg +md=triangle(model(),'../Exp/Square.exp',50000.); +md=setmask(md,'all',''); +md=parameterize(md,'../Par/SquareSheetShelf.par'); +md=setflowequation(md,'SSA','all'); +md.cluster=generic('name',oshostname(),'np',3); + +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); + +% a very special border case that fails the previous version of killberg 20220427 +icebergid = 59; +ids = [141,144,139, 58, 263]; +md.mask.ice_levelset(ids) = 1; +md.mask.ocean_levelset = -md.mask.ice_levelset; +md.mask.ocean_levelset(142) = 1; +md.mask.ocean_levelset(141) = 1; +md.mask.ocean_levelset(143) = 1; +md.mask.ocean_levelset(144) = 1; +md.mask.ocean_levelset(139) = 1; +md.mask.ice_levelset(icebergid) = -1; + +md.levelset.kill_icebergs=1; + +md.timestepping.time_step=10; +md.timestepping.final_time=30; + +%Transient +md.transient.isstressbalance=0; +md.transient.ismasstransport=1; +md.transient.issmb=1; +md.transient.isthermal=0; +md.transient.isgroundingline=0; +md.transient.ismovingfront=1; + +md.calving.calvingrate=zeros(md.mesh.numberofvertices,1); +md.frontalforcings.meltingrate=10000*ones(md.mesh.numberofvertices,1); +md.levelset.spclevelset=md.mask.ice_levelset; + +md=solve(md,'Transient'); + +%Fields and tolerances to track changes +field_names ={'MaskIceLevelsetwithIceberg','MaskIceLevelset1',... + 'MaskIceLevelset2',... + 'MaskIceLevelset3'}; +field_tolerances={1e-11,1e-11,1e-11,1e-11,1e-11,1e-11,1e-11,... + 2e-11,2e-11,2e-11,1e-11,1e-11,1e-11,5e-11,... + 2e-11,2e-11,2e-11,1e-11,1e-11,1e-11,5e-11}; +field_values={... + md.mask.ice_levelset,... + md.results.TransientSolution(1).MaskIceLevelset,... + md.results.TransientSolution(2).MaskIceLevelset,... + md.results.TransientSolution(3).MaskIceLevelset,... + }; +