Index: /issm/trunk-jpl/src/c/cores/hydrology_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 27185) +++ /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 27186) @@ -132,5 +132,5 @@ substep+=1; subtime+=subdt; - averagetime+=subtime*dt; + averagetime+=subtime*subdt; /*Setting substep time as global time*/ femmodel->parameters->SetParam(subtime,TimeEnum); @@ -157,6 +157,6 @@ dtslices=dtslices-remainingslices+(2*remainingslices); subdt=remainingtime/(2*remainingslices); - _printf0_("convergence failed for sub-step "<< substep <<" total number of slice is now "<< dtslices <<" for step "<500){ Index: /issm/trunk-jpl/test/NightlyRun/test906.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test906.m (revision 27186) +++ /issm/trunk-jpl/test/NightlyRun/test906.m (revision 27186) @@ -0,0 +1,92 @@ +%Test Name: SquareNoDynHydrologyDCtimeadapt +md=triangle(model(), '../Exp/Square.exp', 200000.); +md=setmask(md, '', ''); +md=parameterize(md, '../Par/SquareNoDyn.par'); +md.cluster=generic('name', oshostname(), 'np', 1); + +md.transient=deactivateall(md.transient); +md.transient.ishydrology=1; +md.hydrology=(hydrologydc); +md.hydrology=initialize(md.hydrology,md); + +%Hydro Model Parameters +md.hydrology.isefficientlayer=1; +md.hydrology.sedimentlimit_flag=1; +md.hydrology.sedimentlimit=500; +md.hydrology.rel_tol=1.0e-6; +md.hydrology.penalty_lock=10; +md.hydrology.max_iter=100; +md.hydrology.transfer_flag=1; +md.hydrology.unconfined_flag=0; +md.hydrology.leakage_factor=5.0e-10; +md.hydrology.mask_thawed_node=ones(md.mesh.numberofvertices,1); + +%Sediment Parameters +md.hydrology.sediment_porosity=0.4; +md.hydrology.sediment_transmitivity=3.0*ones(md.mesh.numberofvertices,1); +md.hydrology.sediment_thickness=10.; + +%Epl Parameters +md.hydrology.mask_eplactive_node=zeros(md.mesh.numberofvertices,1); +md.hydrology.epl_conductivity=30.; +md.hydrology.epl_initial_thickness=1.0; +md.hydrology.epl_colapse_thickness=1.0e-3; +md.hydrology.epl_thick_comp=0; +md.hydrology.epl_max_thickness=10; +md.hydrology.eplflip_lock=0; + +%Initialisation +md.initialization.sediment_head=md.geometry.base; +md.initialization.epl_head=md.geometry.base; +md.initialization.epl_thickness=ones(md.mesh.numberofvertices,1); + +%Boundary conditions +md.hydrology.spcsediment_head=NaN*ones(md.mesh.numberofvertices,1); +md.hydrology.spcepl_head=NaN*ones(md.mesh.numberofvertices,1); + +%Forcing +stepping=3600./md.constants.yts; +endtime=0.3; +times=[0.:stepping:endtime]; +input=find(abs(times-0.3) > 0.2); +inputval=times*0+100; +inputval(input)=0; + +%basal forcing inputval +md.basalforcings.groundedice_melting_rate=repmat(inputval,md.mesh.numberofvertices,1); +md.basalforcings.groundedice_melting_rate=[[md.basalforcings.groundedice_melting_rate];[times]]; + +%Time +md.timestepping.final_time=endtime; +md.settings.output_frequency=1.0; + +md.timestepping.average_forcing=1; +md.hydrology.step_adapt=1; +md.hydrology.steps_per_step=1; +md.timestepping.time_step=0.1; + + +md=solve(md, 'Transient'); + +field_names={'SedimentWaterHead2', 'EplWaterHead2', 'SedimentHeadResidual2',... + 'SedimentHeadSubstep2', 'EplHeadSubstep2', 'HydrologySubsteps2', 'HydrologySubTime2',... + 'SedimentWaterHead3', 'EplWaterHead3', 'SedimentHeadResidual3'}; +field_tolerances={1e-13, 1e-13, 1e-13,... + 1e-13, 1e-13, 1e-13, 1e-13,... + 1e-13, 1e-13, 1e-13}; +field_values={md.results.TransientSolution(2).SedimentHead,... + md.results.TransientSolution(2).EplHead,... + md.results.TransientSolution(2).SedimentHeadResidual,... + md.results.TransientSolution(2).SedimentHeadSubstep,... + md.results.TransientSolution(2).EplHeadSubstep,... + md.results.TransientSolution(2).HydrologySubsteps,... + md.results.TransientSolution(2).HydrologySubTime,... + md.results.TransientSolution(3).SedimentHead,... + md.results.TransientSolution(3).EplHead,... + md.results.TransientSolution(3).SedimentHeadResidual}; + +for i=1:3 + disp(md.results.TransientSolution(i).EplHead(1)); + disp(md.results.TransientSolution(i).HydrologySubTime); + disp(md.results.TransientSolution(i).HydrologySubsteps); +end Index: /issm/trunk-jpl/test/NightlyRun/test906.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test906.py (revision 27186) +++ /issm/trunk-jpl/test/NightlyRun/test906.py (revision 27186) @@ -0,0 +1,93 @@ +#Test Name: SquareNoDynHydrologyDCtimeadapt +import numpy as np +from model import * +from solve import * +from socket import gethostname +from triangle import triangle +from setmask import setmask +from parameterize import parameterize + +md = triangle(model(), '../Exp/Square.exp', 200000.) +md = setmask(md, '', '') +md = parameterize(md, '../Par/SquareNoDyn.py') + +md.transient = transient.deactivateall(md.transient) +md.transient.ishydrology = True +md.hydrology = hydrologydc() +md.hydrology = md.hydrology.initialize(md) + +#Hydro Model Parameters +md.hydrology.isefficientlayer = 1 +md.hydrology.sedimentlimit_flag = 1 +md.hydrology.sedimentlimit = 500 +md.hydrology.rel_tol = 1.0e-6 +md.hydrology.penalty_lock = 10 +md.hydrology.max_iter = 100 +md.hydrology.transfer_flag = 1 +md.hydrology.unconfined_flag = 0 +md.hydrology.leakage_factor = 5.0e-10 +md.hydrology.mask_thawed_node = np.ones((md.mesh.numberofvertices)) + +#Sediment Parameters +md.hydrology.sediment_porosity = 0.4 +md.hydrology.sediment_transmitivity = 3.0 * np.ones((md.mesh.numberofvertices)) +md.hydrology.sediment_thickness = 10. + +#Epl Parameters +md.hydrology.mask_eplactive_node = np.zeros((md.mesh.numberofvertices)) +md.hydrology.epl_conductivity = 30. +md.hydrology.epl_initial_thickness = 1.0 +md.hydrology.epl_colapse_thickness = 1.0e-3 +md.hydrology.epl_thick_comp = 0 +md.hydrology.epl_max_thickness = 10 +md.hydrology.eplflip_lock = 0 + +#Initialisation +md.initialization.sediment_head = md.geometry.base +md.initialization.epl_head = md.geometry.base +md.initialization.epl_thickness = np.ones((md.mesh.numberofvertices)) + +#Boundary conditions +front = np.where(md.mesh.x == 0)[0] +md.hydrology.spcsediment_head = np.nan * np.ones((md.mesh.numberofvertices)) +md.hydrology.spcepl_head = np.nan * np.ones((md.mesh.numberofvertices)) + +#Forcing +stepping = 3600. / md.constants.yts #hourly forcing +endtime = 0.3 +times = np.arange(0., endtime + stepping, stepping) +inputval = np.where(np.abs(times - 0.3) > 0.2, 0, 100.0) + +#basal forcing inputval +md.basalforcings.groundedice_melting_rate = np.outer(np.ones((md.mesh.numberofvertices)), inputval) +md.basalforcings.groundedice_melting_rate = np.append(md.basalforcings.groundedice_melting_rate, [times], axis=0) + +#Time +md.timestepping.final_time = endtime +md.settings.output_frequency = 1.0 + +md.timestepping.average_forcing = 1 +md.hydrology.step_adapt = 1 +md.hydrology.steps_per_step = 1 +md.timestepping.time_step = 0.1 +md.debug.valgrind = False + +md.cluster = generic('name', gethostname(), 'np', 1) +md = solve(md, 'Transient') + +field_names = ['SedimentWaterHead2', 'EplWaterHead2', 'SedimentHeadResidual2', + 'SedimentHeadSubstep2', 'EplHeadSubstep2', 'HydrologySubsteps2', 'HydrologySubTime2', + 'SedimentWaterHead3', 'EplWaterHead3', 'SedimentHeadResidual3'] +field_tolerances = [1e-13, 1e-13, 1e-13, + 1e-13, 1e-13, 1e-13, 1e-13, + 1e-13, 1e-13, 1e-13] +field_values = [md.results.TransientSolution[1].SedimentHead, + md.results.TransientSolution[1].EplHead, + md.results.TransientSolution[1].SedimentHeadResidual, + md.results.TransientSolution[1].SedimentHeadSubstep, + md.results.TransientSolution[1].EplHeadSubstep, + md.results.TransientSolution[1].HydrologySubsteps, + md.results.TransientSolution[1].HydrologySubTime, + md.results.TransientSolution[2].SedimentHead, + md.results.TransientSolution[2].EplHead, + md.results.TransientSolution[2].SedimentHeadResidual]