![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 27014
- Timestamp:
- 05/19/22 13:30:04 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
r27010 r27014 646 646 int calvinglaw; 647 647 femmodel->parameters->FindParam(&calvinglaw,CalvingLawEnum); 648 IssmDouble mig_max = femmodel->parameters->FindParam(MigrationMaxEnum); 649 IssmDouble dt = femmodel->parameters->FindParam(TimesteppingTimeStepEnum); 650 651 652 /*Get current distance to terminus*/ 653 InputDuplicatex(femmodel,MaskIceLevelsetEnum,DistanceToCalvingfrontEnum); 654 femmodel->DistanceToFieldValue(MaskIceLevelsetEnum,0,DistanceToCalvingfrontEnum); 655 648 656 649 657 if(calvinglaw==CalvingMinthicknessEnum || calvinglaw==CalvingVonmisesEnum){ 650 658 651 659 /*Intermediaries*/ 652 IssmDouble thickness,bed,sealevel ;660 IssmDouble thickness,bed,sealevel,distance; 653 661 IssmDouble min_thickness = femmodel->parameters->FindParam(CalvingMinthicknessEnum); 662 663 /*Loop over all elements of this partition*/ 664 for(Object* & object : femmodel->elements->objects){ 665 Element* element = xDynamicCast<Element*>(object); 666 667 int numnodes = element->GetNumberOfNodes(); 668 Gauss* gauss = element->NewGauss(); 669 Input* H_input = element->GetInput(ThicknessEnum); _assert_(H_input); 670 Input* b_input = element->GetInput(BedEnum); _assert_(b_input); 671 Input* sl_input = element->GetInput(SealevelEnum); _assert_(sl_input); 672 Input* dis_input = element->GetInput(DistanceToCalvingfrontEnum); _assert_(dis_input); 673 674 /*Potentially constrain nodes of this element*/ 675 for(int in=0;in<numnodes;in++){ 676 gauss->GaussNode(element->GetElementType(),in); 677 Node* node=element->GetNode(in); 678 if(!node->IsActive()) continue; 679 680 H_input->GetInputValue(&thickness,gauss); 681 b_input->GetInputValue(&bed,gauss); 682 sl_input->GetInputValue(&sealevel,gauss); 683 dis_input->GetInputValue(&distance,gauss); 684 if(thickness<min_thickness && bed<sealevel && distance<mig_max*dt){ 685 node->ApplyConstraint(0,+1.); 686 } 687 else { 688 /* no ice, set no spc */ 689 node->DofInFSet(0); 690 } 691 } 692 delete gauss; 693 } 694 } 695 else if(calvinglaw==CalvingParameterizationEnum){ 696 697 /*Intermediaries*/ 698 IssmDouble thickness,bed,sealevel,distance; 699 IssmDouble min_thickness = femmodel->parameters->FindParam(CalvingMinthicknessEnum); 654 700 655 701 /*Loop over all elements of this partition*/ … … 662 708 Input* b_input = element->GetInput(BedEnum); _assert_(b_input); 663 709 Input* sl_input = element->GetInput(SealevelEnum); _assert_(sl_input); 710 Input* dis_input = element->GetInput(DistanceToCalvingfrontEnum); _assert_(dis_input); 664 711 665 712 /*Potentially constrain nodes of this element*/ … … 672 719 b_input->GetInputValue(&bed,gauss); 673 720 sl_input->GetInputValue(&sealevel,gauss); 674 if(thickness<min_thickness && bed<sealevel){ 721 dis_input->GetInputValue(&distance,gauss); 722 if(thickness<min_thickness && bed<sealevel && distance<mig_max*dt){ 675 723 node->ApplyConstraint(0,+1.); 676 724 } … … 683 731 } 684 732 } 685 else if(calvinglaw==CalvingParameterizationEnum){ 686 687 /*Intermediaries*/ 688 IssmDouble thickness,bed,sealevel; 689 IssmDouble min_thickness = femmodel->parameters->FindParam(CalvingMinthicknessEnum); 733 else if(calvinglaw==CalvingHabEnum){ 734 735 /*Intermediaries*/ 736 IssmDouble thickness,water_depth,distance,hab_fraction; 690 737 691 738 /*Loop over all elements of this partition*/ … … 693 740 Element* element = xDynamicCast<Element*>(object); 694 741 695 int numnodes = element->GetNumberOfNodes(); 696 Gauss* gauss = element->NewGauss(); 697 Input* H_input = element->GetInput(ThicknessEnum); _assert_(H_input); 698 Input* b_input = element->GetInput(BedEnum); _assert_(b_input); 699 Input* sl_input = element->GetInput(SealevelEnum); _assert_(sl_input); 742 IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum); 743 IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum); 744 745 int numnodes = element->GetNumberOfNodes(); 746 Gauss* gauss = element->NewGauss(); 747 Input* H_input = element->GetInput(ThicknessEnum); _assert_(H_input); 748 Input* bed_input = element->GetInput(BedEnum); _assert_(bed_input); 749 Input* hab_fraction_input = element->GetInput(CalvingHabFractionEnum); _assert_(hab_fraction_input); 750 Input* dis_input = element->GetInput(DistanceToCalvingfrontEnum); _assert_(dis_input); 700 751 701 752 /*Potentially constrain nodes of this element*/ … … 706 757 707 758 H_input->GetInputValue(&thickness,gauss); 708 b_input->GetInputValue(&bed,gauss);709 sl_input->GetInputValue(&sealevel,gauss);710 if(thickness<min_thickness && bed<sealevel){711 node->ApplyConstraint(0,+1.);712 }713 else {714 /* no ice, set no spc */715 node->DofInFSet(0);716 }717 }718 delete gauss;719 }720 }721 else if(calvinglaw==CalvingHabEnum){722 723 /*Intermediaries*/724 IssmDouble thickness,water_depth,levelset,hab_fraction;725 726 /*Get the fraction of the flotation thickness at the terminus*/727 InputDuplicatex(femmodel,MaskIceLevelsetEnum,DistanceToCalvingfrontEnum);728 femmodel->DistanceToFieldValue(MaskIceLevelsetEnum,0,DistanceToCalvingfrontEnum);729 730 /*Loop over all elements of this partition*/731 for(Object* & object : femmodel->elements->objects){732 Element* element = xDynamicCast<Element*>(object);733 734 IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);735 IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);736 737 int numnodes = element->GetNumberOfNodes();738 Gauss* gauss = element->NewGauss();739 Input* H_input = element->GetInput(ThicknessEnum); _assert_(H_input);740 Input* bed_input = element->GetInput(BedEnum); _assert_(bed_input);741 Input* hab_fraction_input = element->GetInput(CalvingHabFractionEnum); _assert_(hab_fraction_input);742 Input* ls_input = element->GetInput(DistanceToCalvingfrontEnum); _assert_(ls_input);743 744 /*Potentially constrain nodes of this element*/745 for(int in=0;in<numnodes;in++){746 gauss->GaussNode(element->GetElementType(),in);747 Node* node=element->GetNode(in);748 if(!node->IsActive()) continue;749 750 H_input->GetInputValue(&thickness,gauss);751 759 bed_input->GetInputValue(&water_depth,gauss); 752 ls_input->GetInputValue(&levelset,gauss);760 dis_input->GetInputValue(&distance,gauss); 753 761 hab_fraction_input->GetInputValue(&hab_fraction,gauss); 754 762 755 if(thickness<((rho_water/rho_ice)*(1+hab_fraction)*-water_depth) && levelset>-300. && levelset<0.){763 if(thickness<((rho_water/rho_ice)*(1+hab_fraction)*-water_depth) && distance<mig_max*dt){ 756 764 node->ApplyConstraint(0,+1.); 757 765 } … … 780 788 781 789 IssmDouble crevasse_threshold = femmodel->parameters->FindParam(CalvingCrevasseThresholdEnum); 782 IssmDouble mig_max = femmodel->parameters->FindParam(MigrationMaxEnum); 783 IssmDouble dt = femmodel->parameters->FindParam(TimesteppingTimeStepEnum); 790 784 791 785 792 for(Object* & object : femmodel->elements->objects){ … … 856 863 surface_input->GetInputValue(&surface,gauss); 857 864 865 /*FIXME: not sure about levelset<0. && levelset>-mig_max*dt! SHould maybe be distance<mig_max*dt*/ 858 866 if((surface_crevasse>surface || crevassedepth>crevasse_threshold*thickness) && bed<0. && levelset<0. && levelset>-mig_max*dt && constraint_nodes[node->Lid()]==0.){ 859 867 local_nflipped++;
Note:
See TracChangeset
for help on using the changeset viewer.
