Question about matrix free method

[abaillod]

Hi,

I am trying to understand what is going on in mtrxhs.f90, and how it is used to evaluate the derivatives of the field w.r.t the geometry degrees of freedom.

@zhisong , would you have a document that explains what is implemented in this part of the code? To me it seems that you evaluate the derivative of the matrix dMA multiplied to the solution array that contains the magnetic field vector potential harmonics. Is this correct?

SPEC/src/mtrxhs.f90

Lines 88 to 136 in 924bcab

	do ii = 1, mn ; mi = im(ii) ; ni = in(ii)
	do ll = 0, lrad
	if (Lcoordinatesingularity) then
	if (ll < mi) cycle ! rule out zero components of Zernike;
	if (mod(ll+mi,2) > 0) cycle ! rule out zero components of Zernike;
	ll1 = (ll - mod(ll, 2)) / 2 ! shrinked dof for Zernike; 02 Jul 19
	else
	ll1 = ll
	end if
	Wte = -two * ni * Tss(ll1,ii) + two * Dzc(ll1,ii)
	Wze = -two * mi * Tss(ll1,ii) - two * Dtc(ll1,ii)
	id = Ate(lvol,0,ii)%i(ll) ; resultA(id) = Wte
	id = Aze(lvol,0,ii)%i(ll) ; resultA(id) = Wze
	if (dBdX%L) cycle ! the derivatives of Lagrange multipliers and helicity w.r.t. interface is zero
	Hte = Ttc(ll1,ii)
	Hze = Tzc(ll1,ii)
	id = Ate(lvol,0,ii)%i(ll) ; resultD(id) = Hte
	id = Aze(lvol,0,ii)%i(ll) ; resultD(id) = Hze
	if( Lcoordinatesingularity .and. ii.eq.1 ) then ; kk = 1
	else ; kk = 0
	endif
	! add Langrange multipliers
	; ; id = Ate(lvol,0,ii)%i(ll) ; resultA(id) = resultA(id)+Lmavalue(lvol,ii) * TTMdata(ll, mi)
	; ; id = Aze(lvol,0,ii)%i(ll) ; resultA(id) = resultA(id)+Lmbvalue(lvol,ii) * TTdata(ll, mi,kk)
	if( ii.gt.1 ) then ; id = Ate(lvol,0,ii)%i(ll) ; resultA(id) = resultA(id)+Lmevalue(lvol,ii) * (- ni * TTdata(ll, mi, 1))
	; ; id = Aze(lvol,0,ii)%i(ll) ; resultA(id) = resultA(id)+Lmevalue(lvol,ii) * (- mi * TTdata(ll, mi, 1))
	else ; id = Ate(lvol,0,ii)%i(ll) ; resultA(id) = resultA(id)+Lmgvalue(lvol,ii) * TTdata(ll, mi, 1)
	; ; id = Aze(lvol,0,ii)%i(ll) ; resultA(id) = resultA(id)+Lmhvalue(lvol,ii) * TTdata(ll, mi, 1)
	endif
	; ; id = Lma(lvol, ii) ; resultA(id) = resultA(id)+Ate(lvol,idx,ii)%s(ll) * TTMdata(ll, mi)
	; ; id = Lmb(lvol, ii) ; resultA(id) = resultA(id)+Aze(lvol,idx,ii)%s(ll) * TTdata(ll, mi,kk)
	if( ii.gt.1 ) then ; id = Lme(lvol, ii) ; resultA(id) = resultA(id)+Ate(lvol,idx,ii)%s(ll) * (- ni * TTdata(ll, mi, 1))
	; ; id = Lme(lvol, ii) ; resultA(id) = resultA(id)+Aze(lvol,idx,ii)%s(ll) * (- mi * TTdata(ll, mi, 1))
	else ; id = Lmg(lvol, ii) ; resultA(id) = resultA(id)+Ate(lvol,idx,ii)%s(ll) * TTdata(ll, mi, 1)
	; ; id = Lmh(lvol, ii) ; resultA(id) = resultA(id)+Aze(lvol,idx,ii)%s(ll) * TTdata(ll, mi, 1)
	endif
	enddo ! end of do ll ;
	enddo ! end of do ii ;

I would like to implement something similar for the new freeboundary version implemented by Stuart - I already implemented the derivatives of the new matrix dMA in the vacuum region, but I don't know how to use these with the matrix free approach.

[zhisong]

Hi,

I am trying to understand what is going on in mtrxhs.f90, and how it is used to evaluate the derivatives of the field w.r.t the geometry degrees of freedom.

@zhisong , would you have a document that explains what is implemented in this part of the code? To me it seems that you evaluate the derivative of the matrix dMA multiplied to the solution array that contains the magnetic field vector potential harmonics. Is this correct?

SPEC/src/mtrxhs.f90

Lines 88 to 136 in 924bcab

	do ii = 1, mn ; mi = im(ii) ; ni = in(ii)
	do ll = 0, lrad
	if (Lcoordinatesingularity) then
	if (ll < mi) cycle ! rule out zero components of Zernike;
	if (mod(ll+mi,2) > 0) cycle ! rule out zero components of Zernike;
	ll1 = (ll - mod(ll, 2)) / 2 ! shrinked dof for Zernike; 02 Jul 19
	else
	ll1 = ll
	end if
	Wte = -two * ni * Tss(ll1,ii) + two * Dzc(ll1,ii)
	Wze = -two * mi * Tss(ll1,ii) - two * Dtc(ll1,ii)
	id = Ate(lvol,0,ii)%i(ll) ; resultA(id) = Wte
	id = Aze(lvol,0,ii)%i(ll) ; resultA(id) = Wze
	if (dBdX%L) cycle ! the derivatives of Lagrange multipliers and helicity w.r.t. interface is zero
	Hte = Ttc(ll1,ii)
	Hze = Tzc(ll1,ii)
	id = Ate(lvol,0,ii)%i(ll) ; resultD(id) = Hte
	id = Aze(lvol,0,ii)%i(ll) ; resultD(id) = Hze
	if( Lcoordinatesingularity .and. ii.eq.1 ) then ; kk = 1
	else ; kk = 0
	endif
	! add Langrange multipliers
	; ; id = Ate(lvol,0,ii)%i(ll) ; resultA(id) = resultA(id)+Lmavalue(lvol,ii) * TTMdata(ll, mi)
	; ; id = Aze(lvol,0,ii)%i(ll) ; resultA(id) = resultA(id)+Lmbvalue(lvol,ii) * TTdata(ll, mi,kk)
	if( ii.gt.1 ) then ; id = Ate(lvol,0,ii)%i(ll) ; resultA(id) = resultA(id)+Lmevalue(lvol,ii) * (- ni * TTdata(ll, mi, 1))
	; ; id = Aze(lvol,0,ii)%i(ll) ; resultA(id) = resultA(id)+Lmevalue(lvol,ii) * (- mi * TTdata(ll, mi, 1))
	else ; id = Ate(lvol,0,ii)%i(ll) ; resultA(id) = resultA(id)+Lmgvalue(lvol,ii) * TTdata(ll, mi, 1)
	; ; id = Aze(lvol,0,ii)%i(ll) ; resultA(id) = resultA(id)+Lmhvalue(lvol,ii) * TTdata(ll, mi, 1)
	endif
	; ; id = Lma(lvol, ii) ; resultA(id) = resultA(id)+Ate(lvol,idx,ii)%s(ll) * TTMdata(ll, mi)
	; ; id = Lmb(lvol, ii) ; resultA(id) = resultA(id)+Aze(lvol,idx,ii)%s(ll) * TTdata(ll, mi,kk)
	if( ii.gt.1 ) then ; id = Lme(lvol, ii) ; resultA(id) = resultA(id)+Ate(lvol,idx,ii)%s(ll) * (- ni * TTdata(ll, mi, 1))
	; ; id = Lme(lvol, ii) ; resultA(id) = resultA(id)+Aze(lvol,idx,ii)%s(ll) * (- mi * TTdata(ll, mi, 1))
	else ; id = Lmg(lvol, ii) ; resultA(id) = resultA(id)+Ate(lvol,idx,ii)%s(ll) * TTdata(ll, mi, 1)
	; ; id = Lmh(lvol, ii) ; resultA(id) = resultA(id)+Aze(lvol,idx,ii)%s(ll) * TTdata(ll, mi, 1)
	endif
	enddo ! end of do ll ;
	enddo ! end of do ii ;

I would like to implement something similar for the new freeboundary version implemented by Stuart - I already implemented the derivatives of the new matrix dMA in the vacuum region, but I don't know how to use these with the matrix free approach.

Yes, it is correct. I should have written some documents for that. But it is basically just computing (15)-(22) of the following page https://princetonuniversity.github.io/SPEC/group__grp__build__matrices.html#ga1c52bb6dbffdb432550aae12b096080f

If you like we can have a face-to-face discussion on Zoom. Currently, I am quite fed up with teaching. The week after next week will be a break so we can do it then.

[abaillod]

Hi Zhisong,

I watched the recording of your SPECtacular talk on December, 10th 2019, and got a general understanding of what you implemented. I still struggle however to see how this can be applied to the new free-boundary approach from Stuart.

I think a zoom meeting between you, @SRHudson, and I would be helpful. Unfortunately, I am on vacation from the 5th till the 23rd of September. It is either this week or we will have to postpone it to the end of September. Let me know what suits you best!