Merge pull request #19 from DJ4Earth/add_ISSM_Archive

Add issm archives

test/runtests.jl

@@ -5,14 +5,46 @@ function searchdir(path,key)
 	filter(x->occursin(key,x), readdir(path))
 end
 
+function compareArchive(id, procedure::Symbol)
+	archive_name = "Archive"*string(id)
+	archive_path = issmdir()*"/test/Archives/"*archive_name*".arch"
+	if procedure===:update
+		# update Archive
+	else
+		# check Archive
+		if isfile(archive_path)
+			for k=1:length(field_names)
+				# Get field and tolerance
+				field=field_values[k];
+				fieldname=field_names[k];
+				tolerance=field_tolerances[k];
+
+				# Compare to archive
+				# Our output is in the correct order (n,1) or (1,1), so we do not need to transpose again
+				archive = archread(archive_path, archive_name*"_field"*string(k))
+				error_diff = (maximum(abs.(archive-field))/(maximum(abs.(archive))+eps(Float64)))
+
+				@test isnan(error_diff) == false
+				@test error_diff < tolerance
+			end
+		else
+			@warn "$archive_name does not exist! Skip the comparison of the results"
+		end
+	end
+end
+
 @time begin
 	@time @testset "Model Struct Tests" begin include("modelstructtests.jl") end
 
 	# test each individual cases, name with test[0-9]*.jl
-	testsolutions = searchdir("./", r"test[0-9]*.jl")
+	testsolutions = searchdir("./", r"test\d+.jl")
 	@time @testset "Model Solution Tests" begin
 		for tf in testsolutions
+			# run the test
 			include(tf)
+			# check the results vs. saved archive
+			testid = match(r"test(\d+).jl", tf).captures[1]
+			compareArchive(testid, :test)
 		end
 	end
 

test/test101.jl

@@ -22,12 +22,12 @@ y     = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","y")
 vx    = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","vx")
 vy    = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","vy")
 index = Int.(archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","index"))
-md.initialization.vx=0 .*InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y,0.0)
-md.initialization.vy=0 .*InterpFromMeshToMesh2d(index,x,y,vy,md.mesh.x,md.mesh.y,0.0)
+md.initialization.vx=InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y,0.0)
+md.initialization.vy=InterpFromMeshToMesh2d(index,x,y,vy,md.mesh.x,md.mesh.y,0.0)
 
-md.materials.rheology_B=1.815730284801701e+08*ones(md.mesh.numberofvertices)
+md.materials.rheology_B=1.815730284801701e08*ones(md.mesh.numberofvertices)
 md.materials.rheology_n=3*ones(md.mesh.numberofelements)
-md.friction.coefficient=20*ones(md.mesh.numberofvertices)
+md.friction.coefficient=0.0*ones(md.mesh.numberofvertices)
 md.friction.p=ones(md.mesh.numberofvertices)
 md.friction.q=ones(md.mesh.numberofvertices)
 
@@ -43,3 +43,10 @@ md.stressbalance.spcvx[pos] .= 0.0
 md.stressbalance.spcvy[pos] .= 0.0
 
 md=solve(md,:Stressbalance)
+
+# Fields and tolerances to track changes
+field_names =["Vx","Vy","Vel"]
+field_tolerances=[4e-13,4e-13,4e-13]
+field_values= [(md.results["StressbalanceSolution"]["Vx"]),
+					(md.results["StressbalanceSolution"]["Vy"]),
+					(md.results["StressbalanceSolution"]["Vel"]) ]