Index: /issm/trunk-jpl/examples/shakti/runme.m =================================================================== --- /issm/trunk-jpl/examples/shakti/runme.m (revision 22371) +++ /issm/trunk-jpl/examples/shakti/runme.m (revision 22372) @@ -3,79 +3,80 @@ if any(steps==1) disp(' Step 1: Mesh'); - - %Generate unstructured mesh on 1,000 m square with typical element edge length of 20 m - md=triangle(model,'./outline.exp',20); - - save MoulinMesh md + + %Generate unstructured mesh on 1,000 m square with typical element edge length of 20 m + md=triangle(model,'./outline.exp',20); + + save MoulinMesh md end if any(steps==2) disp(' Step 2: Parameterization'); - md=loadmodel('MoulinMesh'); - - md=setmask(md,'',''); - - % Run parameterization script to set up geometry, velocity, material properties, etc. - md=parameterize(md,'moulin.par'); - - % HYDROLOGY SPECIFIC PARAMETERIZATION: - % Change hydrology class to Sommers' SHaKTI model - md.hydrology=hydrologysommers(); - - % Define initial water head such that water pressure is 50% of ice overburden pressure - md.hydrology.head = 0.5*md.materials.rho_ice/md.materials.rho_freshwater*md.geometry.thickness + md.geometry.base; - - % Initial subglacial gap height of 0.01m everywhere - md.hydrology.gap_height = 0.01*ones(md.mesh.numberofelements,1); - - % Typical bed bump bump spacing (2m) - md.hydrology.bump_spacing = 2*ones(md.mesh.numberofelements,1); - - % Typical bed bump height (0.1m) - md.hydrology.bump_height = 0.1*ones(md.mesh.numberofelements,1); - - % Define distributed englacial input to the subglacial system (m/yr) - % Change the value 0.0 to add distributed input - md.hydrology.englacial_input = 0.0*ones(md.mesh.numberofvertices,1); - - % Initial Reynolds number (start at Re=1000 everywhere) - md.hydrology.reynolds= 1000*ones(md.mesh.numberofelements,1); - - % Set up atmospheric pressure Type I boundary condition at left edge of - % domain (outflow, i.e. h=zb at x=xmin) - md.hydrology.spchead = NaN(md.mesh.numberofvertices,1); - pos=find(md.mesh.vertexonboundary & md.mesh.x==max(md.mesh.x)); - md.hydrology.spchead(pos)=md.geometry.base(pos); - - save MoulinParam md; + md=loadmodel('MoulinMesh'); + + md=setmask(md,'',''); + + % Run parameterization script to set up geometry, velocity, material properties, etc. + md=parameterize(md,'moulin.par'); + + % HYDROLOGY SPECIFIC PARAMETERIZATION: + % Change hydrology class to Sommers' SHaKTI model + md.hydrology=hydrologysommers(); + + % Define initial water head such that water pressure is 50% of ice overburden pressure + md.hydrology.head = 0.5*md.materials.rho_ice/md.materials.rho_freshwater*md.geometry.thickness + md.geometry.base; + + % Initial subglacial gap height of 0.01m everywhere + md.hydrology.gap_height = 0.01*ones(md.mesh.numberofelements,1); + + % Typical bed bump bump spacing (2m) + md.hydrology.bump_spacing = 2*ones(md.mesh.numberofelements,1); + + % Typical bed bump height (0.1m) + md.hydrology.bump_height = 0.1*ones(md.mesh.numberofelements,1); + + % Define distributed englacial input to the subglacial system (m/yr) + % Change the value 0.0 to add distributed input + md.hydrology.englacial_input = 0.0*ones(md.mesh.numberofvertices,1); + + % Initial Reynolds number (start at Re=1000 everywhere) + md.hydrology.reynolds= 1000*ones(md.mesh.numberofelements,1); + + % Set up atmospheric pressure Type I boundary condition at left edge of + % domain (outflow, i.e. h=zb at x=xmin) + md.hydrology.spchead = NaN(md.mesh.numberofvertices,1); + pos=find(md.mesh.vertexonboundary & md.mesh.x==min(md.mesh.x)); + md.hydrology.spchead(pos)=md.geometry.base(pos); + + save MoulinParam md; end if any(steps==3); disp(' Step 3: Solve!'); md=loadmodel('MoulinParam'); - - md.transient=deactivateall(md.transient); - md.transient.ishydrology=1; - - % Specify that you want to run the model on your current computer - % Change the number of processors according to your machine (here np=4) - md.cluster=generic('np',2); - - % Define the time stepping scheme: run for 90 days with a time step of 1 hr - md.timestepping.time_step=3600/md.constants.yts; % Time step (in years) - md.timestepping.final_time=30/365; - - %Add one moulin with steady input at x=500, y=500 - [a,pos] = min(sqrt((md.mesh.x-500).^2+(md.mesh.y-500).^2)); - time=0:md.timestepping.time_step:md.timestepping.final_time; - md.hydrology.moulin_input = zeros(md.mesh.numberofvertices+1,numel(time)); - md.hydrology.moulin_input(end,:)=time; - md.hydrology.moulin_input(pos,:)=4; - - % Specify no-flux Type 2 boundary conditions on all edges (except - % the Type 1 condition set at the outflow above) - md.hydrology.neumannflux=zeros(md.mesh.numberofelements+1,numel(time)); - md.hydrology.neumannflux(end,:)=time; - - md=solve(md,'Transient'); - save MoulinTransient md + md.transient=deactivateall(md.transient); + md.transient.ishydrology=1; + + % Specify that you want to run the model on your current computer + % Change the number of processors according to your machine (here np=4) + md.cluster=generic('np',2); + + % Define the time stepping scheme: run for 90 days with a time step of 1 hr + md.timestepping.time_step=3600/md.constants.yts; % Time step (in years) + md.timestepping.final_time=30/365; + + %Add one moulin with steady input at x=500, y=500 + [a,pos] = min(sqrt((md.mesh.x-500).^2+(md.mesh.y-500).^2)); + time=0:md.timestepping.time_step:md.timestepping.final_time; + md.hydrology.moulin_input = zeros(md.mesh.numberofvertices+1,numel(time)); + md.hydrology.moulin_input(end,:)=time; + md.hydrology.moulin_input(pos,:)=4; + + % Specify no-flux Type 2 boundary conditions on all edges (except + % the Type 1 condition set at the outflow above) + md.hydrology.neumannflux=zeros(md.mesh.numberofelements+1,numel(time)); + md.hydrology.neumannflux(end,:)=time; + + md.verbose.solution=1; + md=solve(md,'Transient'); + + save MoulinTransient md end Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 22371) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 22372) @@ -287,4 +287,5 @@ ./cores/hydrology_core.cpp\ ./solutionsequences/solutionsequence_hydro_nonlinear.cpp\ + ./solutionsequences/solutionsequence_shakti_nonlinear.cpp\ ./cores/stressbalance_core.cpp\ ./solutionsequences/solutionsequence_stokescoupling_nonlinear.cpp\ Index: /issm/trunk-jpl/src/c/cores/hydrology_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 22371) +++ /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 22372) @@ -113,5 +113,5 @@ femmodel->SetCurrentConfiguration(HydrologySommersAnalysisEnum); InputDuplicatex(femmodel,HydrologyHeadEnum,HydrologyHeadOldEnum); - solutionsequence_nonlinear(femmodel,modify_loads); + solutionsequence_shakti_nonlinear(femmodel); if(VerboseSolution()) _printf0_(" updating gap height\n"); HydrologySommersAnalysis* analysis = new HydrologySommersAnalysis(); Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_shakti_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_shakti_nonlinear.cpp (revision 22372) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_shakti_nonlinear.cpp (revision 22372) @@ -0,0 +1,76 @@ +/*!\file: solutionsequence_nonlinear.cpp + * \brief: core of a non-linear solution, using fixed-point method + */ + +#include "./solutionsequences.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" + +void solutionsequence_shakti_nonlinear(FemModel* femmodel){ + + /*intermediary: */ + Matrix* Kff = NULL; + Matrix* Kfs = NULL; + Vector* ug = NULL; + Vector* uf = NULL; + Vector* old_uf = NULL; + Vector* pf = NULL; + Vector* df = NULL; + Vector* ys = NULL; + + /*parameters:*/ + int max_nonlinear_iterations; + int configuration_type; + IssmDouble eps_res,eps_rel,eps_abs; + + /*Recover parameters: */ + femmodel->parameters->FindParam(&max_nonlinear_iterations,StressbalanceMaxiterEnum); + femmodel->parameters->FindParam(&eps_res,StressbalanceRestolEnum); + femmodel->parameters->FindParam(&eps_rel,StressbalanceReltolEnum); + femmodel->parameters->FindParam(&eps_abs,StressbalanceAbstolEnum); + femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + femmodel->UpdateConstraintsx(); + + int count=0; + bool converged=false; + + /*Start non-linear iteration using input velocity: */ + GetSolutionFromInputsx(&ug,femmodel); + Reducevectorgtofx(&uf, ug, femmodel->nodes,femmodel->parameters); + + //Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate) + InputUpdateFromConstantx(femmodel,converged,ConvergedEnum); + + for(;;){ + /*save pointer to old solution*/ + delete old_uf;old_uf=uf; + delete ug; + + SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel); + CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); + Reduceloadx(pf, Kfs, ys); delete Kfs; + Solverx(&uf, Kff, pf, old_uf, df, femmodel->parameters); + Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete ys; + + convergence(&converged,Kff,pf,uf,old_uf,eps_res,eps_rel,eps_abs); delete Kff; delete pf; delete df; + InputUpdateFromConstantx(femmodel,converged,ConvergedEnum); + InputUpdateFromSolutionx(femmodel,ug); + + /*Increase count: */ + count++; + if(count>=max_nonlinear_iterations){ + _printf0_(" maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded\n"); + converged = true; + } + if(converged) break; + } + + if(VerboseConvergence()) _printf0_("\n total number of iterations: " << count << "\n"); + + /*clean-up*/ + delete uf; + delete ug; + delete old_uf; +} Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h (revision 22371) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h (revision 22372) @@ -14,4 +14,5 @@ void solutionsequence_thermal_nonlinear(FemModel* femmodel); void solutionsequence_hydro_nonlinear(FemModel* femmodel); +void solutionsequence_shakti_nonlinear(FemModel* femmodel); void solutionsequence_nonlinear(FemModel* femmodel,bool conserve_loads); void solutionsequence_newton(FemModel* femmodel);