Check resolution of BCs in adaptive loop

ralfHielscher · Nov 20, 2017 · 9a75c18 · 9a75c18
1 parent 6d09efd
commit 9a75c18
Show file tree

Hide file tree

Showing 3 changed files with 26 additions and 14 deletions.
diff --git a/@chebop2/solvepde.m b/@chebop2/solvepde.m
@@ -100,22 +100,27 @@
 Resolved_y = 0; 
 Resolved = Resolved_x & Resolved_y;
 
+% Should we check the resolution of the BCs?
+bctype = chebop2.checkBC(N, 2, 2);
+
 while ( ( ~Resolved ) && ( m < maxDiscretise_y ) &&...
                          ( n < maxDiscretise_x ) )
 
     % Solve PDE, return an m x n matrix of coefficients:
-    X = chebop2.denseSolve(N, f, n, m);
+    X = chebop2.denseSolve(N, f, m, n);
+    old_m = m;
+    old_n = n;
 
     if ( adaptive_y ) 
         % Resolved in y-direction?
-        [Resolved_y, m] = resolveCheck(X, m, tol);
+        [Resolved_y, m] = resolveCheck(X.', old_m, old_n, tol, N.lbc, N.rbc, bctype);
     else
         Resolved_y = 1; 
     end
 
     if ( adaptive_x ) 
         % Resolved in x-direction?
-        [Resolved_x, n] = resolveCheck(X.', n, tol);
+        [Resolved_x, n] = resolveCheck(X, old_n, old_m, tol, N.dbc, N.ubc, bctype);
     else 
         Resolved_x = 1; 
     end
@@ -147,17 +152,24 @@
 
 end
 
-function [Resolved, newDisc] = resolveCheck(coeffs, oldDisc, tol)
+function [Resolved, newDisc] = resolveCheck(coeffs, m, n, tol, lbc, rbc, bctype)
 % Basic Resolution check:
 
 tail = max(abs(coeffs(:,end-8:end)));
-Resolved = all(tail < 20*oldDisc*tol);
+Resolved = all(tail < 20*m*tol);
+
+% Check resolution of the Dirichlet BCs
+if ( bctype == 1 )
+    lval = (-1).^(0:n-1) * coeffs; lbc_cfs = chebcoeffs(lbc, m);
+    rval = ones(1,n) * coeffs;     rbc_cfs = chebcoeffs(rbc, m);
+    Resolved = Resolved & norm(lval' - lbc_cfs) < tol & norm(rval' - rbc_cfs) < tol;
+end
 
 % If unresolved, then increase the grid:
 if ( ~Resolved )
-    newDisc = 2^(floor(log2(oldDisc))+1) + 1;
+    newDisc = 2^(floor(log2(m))+1) + 1;
 else
-    newDisc = oldDisc;
+    newDisc = m;
 end
 
 end
@@ -167,7 +179,7 @@
 
 % Increase tolerance on large grids:
 if ( max(m, n) > 250 )
-    tol = max(tol, 1e-10);
+    tol = max(tol, 1e-11);
 end
 
 end

diff --git a/tests/chebop2/test_adaptivity.m b/tests/chebop2/test_adaptivity.m
@@ -16,27 +16,27 @@
 N.rbc = @(t,u) [u ; diff(u)];
 N.lbc = 0;
 u1 = mldivide(N, 0, 200, inf);
-[m, n] = length( u1 ); 
+[n, m] = length( u1 );
 pass(1) = ( m == 200 ); 
 
 u2 = mldivide(N, 0, 200, 200);
-[m, n] = length( u2 ); 
+[m, n] = length( u2 );
 pass(2) = ( m == 200 ); 
 pass(3) = ( n == 200 ); 
 
 u3 = mldivide(N, 0, inf, 200);
-[m, n] = length( u2 ); 
+[n, m] = length( u2 );
 pass(4) = ( n == 200 ); 
 
 % check small values of n: 
 u4 = mldivide(N, 0, 10, 10);
-[m, n] = length( u4 ); 
+[n, m] = length( u4 );
 pass(5) = ( m == 10 ); 
 pass(6) = ( n == 10 );  
 
 % check small values of n: 
 u4 = mldivide(N, 0, 6, 10);
-[m, n] = length( u4 ); 
+[n, m] = length( u4 );
 pass(7) = ( n == 10 ); 
 pass(8) = ( m == 6 );
 

diff --git a/tests/chebop2/test_domain.m b/tests/chebop2/test_domain.m
@@ -5,7 +5,7 @@
 if ( nargin < 1 ) 
     prefs = chebfunpref(); 
 end 
-tol = 100*prefs.techPrefs.chebfuneps; 
+tol = 1e3*prefs.techPrefs.chebfuneps;
 
 d = [-2 2 -2 2];
 N = chebop2(@(u) diff(u, 2, 1) + diff(u, 2, 2), d);