New regularization over the current potential

[zhucaoxiang]

@landreman @ejpaul As both of you know, I am working on using REGCOIL to help design permanent magnets. I implemented a new regularization term over the integral of current potential squared. Details are shown in the following figure.

I have pushed a new branch phi_regularize and the key part is as following.

regcoil/regcoil_build_matrices.f90

Lines 602 to 631 in 17b94d4

	case (regularization_term_option_chi2_Phi)
	! ------------------------------------------------------------------
	! current potential matrix and RHS:
	call system_clock(tic)
	! Here we carry out matrix = dtheta*dzeta*(basis ^ T) * basis
	! A = basis_functions
	! B = basis_functions
	! C = matrix
	M = num_basis_functions ! # rows of A^T
	N = num_basis_functions ! # cols of B
	K = ntheta_coil*nzeta_coil ! Common dimension of A^T and B
	LDA = K ! Would be M if not taking the transpose.
	LDB = K
	LDC = M
	TRANSA = 'T' ! DO take a transpose!
	TRANSB = 'N'
	BLAS_ALPHA = dtheta_coil*dzeta_coil
	BLAS_BETA=0
	call DGEMM(TRANSA,TRANSB,M,N,K,BLAS_ALPHA,basis_functions,LDA,basis_functions_sqrt_N_coil,LDB,BLAS_BETA,matrix_regularization,LDC)
	call system_clock(toc)
	if (verbose) print *,"matmul for matrix_regularization:",real(toc-tic)/countrate,"sec."
	call system_clock(tic)
	RHS_regularization = matmul(reshape(secular_current_potential, (/ ntheta_coil*nzeta_coil /)), basis_functions_sqrt_N_coil)
	call system_clock(toc)
	if (verbose) print *,"Matmul for RHS_regularization:",real(toc-tic)/countrate,"sec."

---

Initial results

I have a test on the new current potential regularization, using a NCSX saddle coils only case.

Regularization on K

I plotted out the chi2_Phi from the output regularized with the current density. The current potential is actually having the same pattern, although it is not directly penalized.

Here are plots of Phi and K.

Regularization on Phi

Here are the results of Phi regularization.

Discussions

The L-curves look similar (and sensible), although regularization on Phi seems to get more localized distribution, both for Phi and K. This is not quite as I expected.
Since you two are definitely more familiar with the source codes, could you please help me check if I was implementing the correct regularization? If yes, I can create a pull-request and merge the commits.

[zhucaoxiang]

I will close this because it appears the chi^2_K is outperforming chi^2_Phi. Because the chi^2_K is regularizing the gradient of Phi, it tends to provide more smooth solutions.