Commit d24427fb by Lars Stietz

### sweep_Mat function implemented

parent a750b6df
 ... ... @@ -48,46 +48,54 @@ def asmSlowMat(dt): Q = Q + np.diag(dt[i]*np.ones(sz-i),-i) return Q # function for SDC sweep in one iteration step def sweep(u,d_tau,Smat,Qmat,f): """ adjustment for strt/nd=False still needed """ num_nodes = len(u) u_old = deepcopy(u) for m in range(1,num_nodes-1): u[m] = 1/(1-f(1)*d_tau[m])*(u[m-1]-d_tau[m]*f(u_old[m]) + np.dot(Smat[m+1,1:],f(u_old))) u[num_nodes] = u[0] + np.dot(Qmat[-1,1:],u) u[num_nodes-1] = u[0] + np.dot(Qmat[-1,1:],f(u)) return u # function for SDC sweep with matrices def sweepMat(u,d_tau,Qmat,Qs,f): u_old = deepcopy(u) u = np.linalg.solve(np.eye(np.shape(Qs)[0])-f(1)*Qs,u[0]-f(1)*Qs@u_old+f(1)*Qmat[1:,1:]@u_old) u[-1] = u[0] + np.dot(Qmat[-1,1:],f(u)) return u # function to calculate the residual def res(u, Qmat,u_0): """ adjustment for strt/nd=False still needed """ r = u_0 + Qmat[1:,1:]@f(u)-u return r solving = 1 #1 solves straight forward, 1 looks at accuracy, 2 compares matrix approach to straight forward solving = 1 # 0<->accuracy, 1<->solving, 2<->matrix sweep comparison if solving == 1: # SDC straight forward approach """ res = np.zeros(num_nodes) k_max = 20 res_max = np.amax(u[0] +coll.Qmat[2 :,:]@u) # SDC straight forward approach k_max = 100 res_2 = 2 res_it = np.empty(0) j = 0 error = 1 error_it = np.empty(0) """ for it in range(250): #while res_max>1e-15 or error>1e-15: sweep(u,coll.delta_m,coll.Smat,coll.Qmat,f) Q_slow = asmSlowMat(coll.delta_m) #while res_2 > 1e-14: for it in range(k_max): #u = sweep(u,coll.delta_m,coll.Smat,coll.Qmat,f) u = sweepMat(u,coll.delta_m,coll.Qmat,Q_slow,f) r = res(u.real,coll.Qmat,u_0) """ for i in range(num_nodes-1): res[i] = u[0].real + np.dot(coll.Qmat[i][1:-1],f(u[1:-1]).real)-u[i].real res_max = np.amax(abs(res)) error = np.linalg.norm(u_old.real-u.real) error_it = np.append(error_it,error) res_it = np.append(res_it,res_max) """ res_2 = np.linalg.norm(r) res_it = np.append(res_it,res_2) plt.figure(1) ... ... @@ -100,8 +108,8 @@ if solving == 1: plt.legend() plt.figure(2) plt.title('Update ans Residual') plt.semilogy(error_it,label= 'update difference (2-Norm)') plt.semilogy(r,label='residual') #plt.semilogy(error_it,label= 'update difference (2-Norm)') plt.semilogy(res_it,label='residual') plt.xlabel('iteration step') plt.legend() ... ... @@ -111,9 +119,7 @@ elif solving == 0: kons = np.zeros(it) for j in range(it): u_old = deepcopy(u) for i in range(num_nodes): u[i+1] = 1/(1-l*coll.delta_m[i])*(u[i]-coll.delta_m[i]*f(u_old[i+1])+np.dot(coll.Smat[i+1][1:],f(u_old[1:]))) u = sweep(u,coll.delta_m,coll.Smat,coll.Qmat,f) kons[j] = np.log(abs((y(coll.nodes[-1])-u[-1])))/np.log(np.amax(coll.delta_m)) if j>0: print(kons[j]-kons[j-1]) ... ... @@ -151,11 +157,3 @@ elif solving == 2: plt.legend() plt.show() \ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!